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ABSTRACT 


As  an  alternative  to  the  LU  matrix  factorization,  we  consider  a  factoriza¬ 
tion  that  uses  the  lower  triangular  part  of  the  original  matrix  as  one  factor  and 
computes  the  other  factors  as  a  product  of  rank-one  update  matrices. 

Under  some  non-singularity  assumptions,  an  m  xm  matrix  A  can  be  factorized 
as  EmEm-i  •  •  •  E2  Ai  where  Ai  is  the  lower  triangular  part  of  A  and  Et  is  a  rank- 
one  update  matrix  of  the  form  I  -j-  with  vjt  a  column  vector  and  u;*  a  row 

vector.  The  vector  is  the  column  of  A  — Ai.  If  =  0,  then  E^  =  I  may  be 
omitted  from  the  factorization.  Otherwise,  the  row  vector  uft  must  be  computed. 

After  reviewing  and  improving  the  time  complexity,  the  requirements,  the 
stability  and  the  efficiency  of  this  method,  we  derive  a  stable  factorization  algo¬ 
rithm  which  we  implement  in  FORTRAN77  within  the  framework  of  the  simplex 
algorithm  for  linear  programming. 

A  comparison  of  our  numerical  results  with  those  obtained  through  the  code 
MINOS  5.3  indicate  that  our  method  may  be  more  efficient  than  an  ordinary 
LU  decomposition  for  some  matrices  whose  order  ranges  between  28  and  1481, 
especially  when  these  matrices  are  almost  triangular. 
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INTRODUCTION 


The  most  widely  used  matrix  factorization,  the  LU  factorization,  amounts 
to  the  computation  of  two  triangular  factors,  one  of  which  can  be  regarded  as 
a  product  of  rank-one  update  matrices.  As  an  alternative,  George  B.  Dantzig 
(1985)  has  proposed  another  factorization  that  uses  the  lower  triangular  part  of 
the  original  matrix  as  one  factor  and  computes  the  other  factor  as  a  product  of 
rank-one  update  matrices. 

Under  some  non-singularity  assumptions,  anmxm  matrix  A  can  be  factorized 
as  ■  •  •  E2A1  where  Ai  is  the  lower  triangular  part  i^f  A  and  Ejt  is  a  rank- 

one  update  matrix  of  the  form  I  -f  VkUJk  with  v*  a  column  vector  and  uJk  a  row 
vector. 

The  vector  vjt  is  the  A:**  column  of  A  —  Aj.  If  vt  =  0,  then  Efc  =  I  may  be 
omitted  from  the  factorization.  Otherwise,  the  row  vector  UJk  defining  E^  can  be 
obtained  by  solving  a>jtEjt_i  . . .  E2  Ai  =  where  is  the  unit  row  vector. 

Once  the  auxiliary  vectors  have  been  computed,  any  system  Ax  =  b  or 
ttA  =  7  can  be  solved  through  one  sparse  triangular  system  involving  Ai  and 
s  rank-one  updates  involving  the  matrices  Efc  for  which  column  fc  is  a  spike  (i.e. 
V*:  ^0). 

However,  that  factorization  may  break  down,  for  instance  on  a  matrix  whose 
diagonal  contains  a  zero  element.  In  addition,  even  if  the  factorization  exists,  it 
will  often  be  unstable.  In  this  thesis,  we  show  how  to  overcome  the  problems  of 
existence  and  stability.  We  present  a  factorization  method  theoretically  applicable 
to  any  non-singular  matrix.  The  numerical  results  presented  in  the  last  chapter 
indicate  that  our  method  may  be  more  efficient  than  an  ordinary  LU  factorization 
for  some  matrices  whose  order  ranges  from  28  to  1481,  and  whose  number  of 
elements  ranges  from  72  to  6344. 

In  Chapter  1,  we  indicate  how  to  compute  the  vectors  recursively,  and 
how  to  solve  systems  such  as  Ax  =  b  and  ttA  =  7  using  these  vectors,  assuming 
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the  non-singularity  of  A  and  of  some  of  its  submatrices. 

In  Chapter  2,  we  streamline  the  method  described  in  Chapter  1  and  reduce 
its  computational  complexity  to  the  same  level  as  that  of  the  LU  factorization. 

In  Chapter  3,  we  study  the  existence  and  the  stability  of  the  factorization. 
We  introduce  the  vectors  fTk  (2  <  k  <  m),  a  normalized  form  of  the  vectors 
tJfc  {2  <  k  <  m).  K  A  is  non-singular,  a  permutation  Q  of  the  columns  of  A 
makes  AQ  and  its  leading  principal  submatrices  non-singular  and,  in  practice, 
well  conditioned.  We  indicate  how  to  compute  simultaneously  the  permutation  Q 
and  the  vectors  representing  the  factorization  of  AQ. 

In  Chapter  4,  we  present  some  algorithms  inspired  by  Hellerman  and  Rarick 
(1971),  that  permute  the  rows  and  columns  of  A  in  order  to  reduce  the  number 
of  spikes  (i.e.  the  columns  that  have  nonzero  entries  above  the  diagonal).  The 
method  is  more  efficient  if  there  are  fewer  spikes  and  if  the  spikes  are  shifted 
towards  the  left  because  ufk  has  at  most  k  nonzero  components. 

In  Chapter  5,  we  propose  an  algorithm  that  reorders  the  rows  and  columns  of 
A  while  computing  the  vectors  cr^.,  so  as  to  reduce  the  number  of  spikes  without 
compromising  the  numerical  stability  of  the  factorization. 

In  Chapter  6,  we  describe  a  FORTRAN  implementation  of  the  algorithm 
within  the  framework  of  linear  programming.  We  use  our  own  set  of  factorization 
routines  in  the  optimization  code  MINOS  5.3  of  Murtagh  and  Saunders  (1987)  on 
a  set  of  51  test  problems  from  netlib  (Gay,  1985).  With  the  options  chosen,  our 
method  achieves  faster  running  times  than  the  original  MINOS  code  on  about  a 
third  of  the  problems.  Unfortunately,  it  also  takes  up  to  four  times  as  long  as 
MINOS  on  some  problems. 
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CHAPTER  1.  THE  BASIC  METHOD 


1.1  Introduction 

Under  sonae  non- singularity  assumptions,  anmxm  matrix  A  can  be  factorized 
into  EmEm  — 1  •  •  •  E2A1  where  Ai  is  the  lower  triangular  part  of  A  amd  Ejt  is  a 
rank-one  update  matrix  of  the  form  I  -I-  VfcW/t  with  Vjt  a  column  vector  and  uJk  a 
row  vector. 

The  vector  V/t  is  the  column  of  A  —  Ai.  If  Vfc  =  0,  then  Et  =  I  may  be 
omitted  from  the  factorization.  Otherwise,  the  row  vector  must  be  computed. 

In  this  Chapter,  we  describe  the  factorization  and  provide  a  natural  method 
to  compute  the  vectors  u)k  (2  <  k  <  m)  and  to  solve  systems  such  as  Ax  =  b  or 
ttA  =  y  using  this  factorization. 


1.2  Definitions  and  Notation 

Let  m  be  a  positive  integer.  Let  R  denote  the  real  numbers.  Unless  otherwise 
specified,  matrix  denotes  an  element  of  R”*^"*,  column  vector  denotes  an  element 
of  R"*^i  and  row  vector  denotes  an  element  of  R^^"*.  Row  vectors  are  usually 
represented  by  Greek  letters. 

A  given  matrix  A  can  be  decomposed  into  its  lower  triangiileu:  part  Aj  (in¬ 
cluding  the  diagonal)  and  its  strictly  upper  tri^mgula^  part,  which  can  in  turn  be 
broken  down  according  to  its  entries  in  columns  2, . . . ,  m  into  m  —  1  matrices  of 
the  form  for  2  <  A:  <  m,  where  u*  is  the  imit  vector.  Then 

A  =  Ai  -h  -I-  vsuf  -I-  ••  + 

For  2  <  A:  <  m,  let 

.A  -  'T 

A*  =  Afc_i  +  v*ujj 

=  Ai  -I-  V2u|’  -!-•••  +  VfcU^. 

If  Vfc  ^0,  column  k  is  called  a  spike  of  A. 
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In  this  chapter,  we  assume  that  the  matrices  A*  {I  <  k  <  m)  are  non¬ 
singular  (cf.  Chapter  3).  We  can  then  turn  the  above  decomposition  of  A  into  a 
factorization.  For  2  <  k  <  m,  let 

I  ^ 

I  Efc  =  I  +  VfcWt 

Then 

A*  =  A*_i  +  v*uj’ 

=  (I  +  VfcU^A^iJ  At_i 

=  (I  +  VfcU>*)  Afc_i 

=  EfcAfc_i 

=  E*Ejk-i .  • .  E2Ai 

and  in  particular,  A  =  En,Em-i  •  •  ■  EzAi.  Each  matrix  Ej^  is  a  rank-one  update 

matrix.  Note  that  we  could  also  have  factored  Ai  on  the  left  and  obtained  a 

factorization  of  the  form  A1F2  . . .  Fjk-iF*  where  each  matrix  F^  is  also  a  rank- 
one  update  matrix. 
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1.3  Solution  of  Rank-One  Update  Systems 


Let  E  =  I  VLj  and  ^  =  1  4-  u>v.  The  following  Lemmas  are  well  known 
properties  of  rank-one  update  matrices  (Golub  &  Van  Loan,  1983). 

Lemma  1.3.1 

(  d  =  det  E 

1  E-'  =  I  - 

Lemma  1.3.2  The  solutions  to  the  elementzury  linear  systems  Ex  =  b  and 
ttE  =  7  are  given  by  the  following  rank-one  update  formulas: 

^  =  1  -t-  u>v 

<  X  =  b  -  9~^  (tt^b)  V 

?r  =  7  -  (7v)  u} 

Once  LJ  zmd  6  have  been  computed,  the  number  of  multiplications  or 
divisions  required  in  either  update  is  at  most  equal  to  the  number  of  nonzero 
components  in  v,  plus  the  number  of  nonzero  components  in  plus  one  (corre¬ 
sponding  to  the  factor  By  construction,  the  vectors  v*  and  defined  in 

Section  1.2  have  at  most  —  1  and  k  nonzero  components  respectively.  Therefore, 
solving  the  system  EfcX  =  b  or  jrE/t  =  7  requires  at  most  2k  multiplica¬ 
tions  or  divisions.  The  same  upper  bound  holds  for  the  number  of  additions  or 
subtractions. 


1.4  Computation  of  the  Factors  E*  (2  <  A:  <  m) 

For  our  purposes,  knowing  the  factor  E*  =  I  -f-  through  ohe  vec¬ 

tors  Vfc  amd  u>k  is  sufficient.  Since  the  vectors  vt  are  given,  we  need  only 
compute  the  auxilieiry  vectors  cofc  {2  <  k  <  m)  defined  by  =  uj,  i.e. 
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UfcEfc-i  . . .  E2  Ai  =  uj .  This  can  be  done  recursively  by  solving  the  systems 

UJ2A1  = 
a>3E2Ai  = 

< 

ttJfcEfc-i . . .  E2  Ai  = 

^  — lEnj  — 2  •  •  •  E2 Aj  = 

in  that  order.  To  solve  the  system  OJtEjt-i . . .  E2A1  =  u^,  we  can  solve  the 
triangular  system  ffiAj  =  by  backward  substitution  and  then  successively 
solve  k  —  2  rank -one  update  systems  of  the  form  ir/E/  =  iti-i  for  2  <  /  <  it  —  1. 
By  Lemma  1.3.2,  these  systems  are  equivalent  to 

jr,  =  *  (Jr/_iv/)  for  2  <  1  <  k  ~  1. 

Once  oJk  =  IT k-i  is  known,  we  can  obtain  9k  by 

=  1  -I-  a;*Vfc. 

Since  the  last  m  —  k  components  of  iTi  and  its  successive  rank-one  updates 
(including  are  zero,  all  these  systems  are  of  dimension  at  most  k  for  practical 
purposes.  The  maiximum  number  of  multiplications  required  to  compute  (u;*,  6k) 
by  this  method  is 

*  — 1  o 

ijfc(jt-|.l)  +  Y^21  -f  Jb-1  w  |jfc2 +fc 

1=2 

and  the  maximum  number  of  multipUcations  required  to  factorize  A  is 

k=2 

Example: 

Consider  the  3x3  matrix  A  =  I  3  5  7  I  introduced  in  Section  1.2. 

\4  9  2j 
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The  vector  u>2  characterizing  E2  is  given  by  u;2Ai  =  u^. 


UJ2 


/8  . 

3  5 
\4  9 


(0  1  0) 


^^2  =  (-^  i  *) 

^2  =  1  +  (-:^  5 


37 

40 


The  vector  characterizing  E3  is  given  by  u>>3E2Ai  =  u^. 

=(001) 

tJ3E2=(--^  i) 


/  i  —  3  9  I  \  40 

^*^3  —  40  10  2  ^  37 


f-J.  -±  1) 

t  40  10  2  / 


/  .  —  7  34  1  \ 

^3  —  (  74  37  2  ) 


(- 


40 


T  * 


^3  =  1  +  (^  -U  i)  7  =  - 


—  180 
37 


1.5  Solution  of  Ax  =  b 

To  solve  the  system  Ax  =  b  given  the  non-singular  factorization  A 
Em  •  •  •  E2A1,  the  natural  method  is  to  solve  successively  the  systems 

f  Embm  — 1  =  bm 


Etbfc_i  =  bfc 

E2bi  =  b2 

AjX  =  bj 
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for  bOT-i,...,bi  and  x,  having  let  b^  =  b.  Applying  Lemma  1.3.2  to  the 
non-singular  matrices  Efc  (m  >  *  >  2)  ,  we  can  solve  the  first  m  -  1  systems 
according  to  the  ramk-one  update  formulas 

b*-i  =  bfc  -  6^^  (ujfcbfc)  V*  for  m>k>2 

aind  then  solve  the  triangular  system  AiX  =  bi  by  forward  substitution. 

The  maximum  number  of  multiplications  required  by  this  method  is 
2 

2k  -I-  ^m(m  -t- 1)  «  |m^  +  |m. 

k—m 


Example: 

/8  1  6\ 

Consider  the  3x3  system  Ax  =  b  where  A  =  I  3  5  7  I  and  b  = 

\4  9  2/ 

/24\ 

I  14  j .  The  system  can  be  written  E3E2A1X  =s  b.  The  vectors  ba,  ba,  bj  and 

\-sJ 

X  are  the  following; 


/24\ 
ba  =  I  14 
\-8/ 
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1.6  Solution  of  ttA  =  7 

To  solve  the  system  irA  =  7  given  the  non-singulax  factorization  A  = 
Em  •  ■  ■  E2A1,  the  natural  method  is  to  solve  successively  the  systems 

7riAi=7 
7r2E2  =  ■JTi 

TTfcEi  =  Kk-i 

'TTmErn  ~  ^m  — 1 

for  JTi ,  ir2, . . . ,  TTm  =  JT-  The  triangular  system  itiAj  =7  can  be  solved  by 
backward  substitution.  Then,  applying  Lemma  1.3.2  to  the  non-singular  matrices 
Efc  ,  we  can  solve  the  last  m  —  1  systems  according  to  the  rank-one  update  formulas 

Kk  =  iTfc-i  -  ^  (-TTk-iyk)  u^it  for  2<k<  m. 

The  maximum  number  of  multiplications  required  by  this  method  is 

m 

im(m -I- 1)  +  -f- |m. 

fc=2 


Example: 

8  1  6\ 

3  5  7  1  and  7  = 

4  9  2/ 

(6  0  —6).  The  system  can  be  written  2rE3E2Ai  =  7.  The  vectors  Wi,  ^2,  Ws 

zmd  TT  are  the  following: 


Consider  the  3x3  system  wA  =  7  where  A  = 


(6  0  -6) 


’Ti  =  (^  ¥  -3) 
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^3 


V  40 

f  -3) 

40  / 

37  V 

9  27 

40  5 

(■^ 

V  37 

198  0  \ 

37  ^  ) 

(  -2- 
V  37 

198  0  \ 

TT  ) 

^  180 

(  9  198 

V  37  37 

(1 

-2  1) 

i-^0  I  *) 


(  1.  -li  i) 

V  74  37  2  / 


2r  =  (1  -2  1) 


1.7  Fundamental  Observation 

If  column  I  is  not  a  spike  (i.e,  if  v/  =  0),  then  E(  =  1  may  be  omitted 
from  the  factorization  of  A  and  the  rank-one  update  corresponding  to  E;  may 
be  skipped  when  solving  Ax  =  b,  irA  =  7  or  even  u>fcAj^-i  =  u^.  This 
simplification  significantly  reduces  the  size  of  the  computations  when  the  matrix 
A  is  Izirge  and  has  few  spikes. 


Example: 

This  example  shows  how  to  factorize  a  5  x  5  matrix  whose  columns  2  and  4 
are  not  spikes,  and  how  to  solve  systems  like  Ax  =  b  or  irA  =  'y. 


A  = 


/  1 

0 

1 

0 

1\ 

0 

0 

1 

0 

0 

1 

-1 

0 

1 

0 

1 

b  = 

6 

0 

1 

1 

1 

0 

6 

\  1 

1 

0 

1 

U/ 

7  =  (1  3  3  2  3) 
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Ai  = 


-1 
0 
V  1 


A2  =  Ai 

/  1 
0 

As  =  —1 

V  1 

A4  =  A3 

/  1 
0 

-1 


As  = 


0 
\  1 


*  *  *  *  \ 

1  *  *  * 

0  1** 
111* 
10  11/ 


0  1  *  *\ 
10** 

0  1** 
111* 
10  11/ 


0  1  0  i\ 
10  0  0 
0  10  1 
1110 
10  11/ 


Vo  =  0 


Eo 


V3  = 


/1\ 

0 

v:/ 


V4  =  0 


Vs  = 


1 

0 

\*/ 


E4 


•  The  vector  W3  characterizing  E3  is  given  by  u^sAi  =  Uj  . 


U>3 


1 

0 

-1 

0 

,  1 


*  *  *  *\ 

1  *  *  * 

0  1** 
111* 
10  11/ 


=  (0  0  1  0  0) 


U>3  =  (  1  0  1  *  *  ) 

^3  =  1  +  (1  0  1  * 


2 


T 

•  The  vector  u>5  characterizing  E5  is  given  by  uJsEsAi  =  U5  . 


W5E3 


/  1 
0 

-1 

0 

V  1 


*  ♦  ♦  *\ 

1  ♦  ♦  * 

0  1** 
111* 
10  11/ 


=  (0  0  0  0  1) 
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U>5E3  =  (0  0  1  -1  1) 


=  (0  0  1  -1  1)  -  i  (0  0  1  -1  1)  *  (1  0  1  ♦  *) 


a;5  =  (0  0  1  -1  1) 


^5  =  1  +  (0  0  1  -1  1)  1  =2 

0 


•  The  system  Ax  =  b  can  be  written  E5E3A1X  =  b.  The  vectors  bj,  bs,  bi 
and  X  are  the  following: 


bs  = 


b3  = 


i  (0  0  1  -1  1)  6 

6 


bi  = 


i  (1  0  1  ♦  ♦)  2 

6 


1  *  ♦  ♦ 

0  1***1 
-1  0  1  *  ♦  I  X 
0  111*1 
1  10  11/ 
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•  The  system  ttA  =  7  can  be  written  flrEsEsAi  =  7.  The  vectors  Wi,  Vj, 
TTs  and  IT  are  the  following: 


^1 


/I  *  *  *  ♦x 

0  1  *  *  *  ' 

-1  0  1  ♦  * 

0  1  1  1  * 

\  1  1011/ 


=  (1  3  3  2  3) 


TTi  =  (2  1  4  -1  3) 


ITS  =  (2  1  4  -1  3) 
=(113-13) 
TTs  =  (1  1  3  -1  3) 
=(11111) 


(2  14-1  3) 


♦ 

♦ 

{*/ 


(10  1**) 


0 

(113-1  3) 

1 

0 

[*J 

(0  0  1-1  1) 


TT  =  (1  1  1  1  1) 
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CHAPTER  2.  THE  STREAMLINED  METHOD 


2.1  Introduction 

The  basic  method  of  Chapter  1  is  not  the  most  efficient  way  to  compute 
and  use  the  factorization  A  =  E^Em-i  •••EaAi.  In  this  Chapter,  we  present 
a  streamlined  method  that  yields  the  same  results  as  the  basic  method  without 
explicitly  solving  the  intermediate  system  involving  the  triangular  matrix  Ai .  This 
decreases  by  about  33%  the  upper  bound  on  the  number  of  multiplications  required 
to  compute  the  factorization  of  A,  the  solution  of  Ax  =  b  or  the  solution  of 
irA  =  7. 


2.2  Solution  of  Ax  =  b 

The  basic  method  explained  in  Section  1.5  can  be  streamlined  by  observing 
that  Xk,  the  component  of  x,  is  computed  when  the  rank-one  update  system 
Efcbik-i  =  bfc  is  solved.  Consider  the  following  Lemmas: 

Lemma  2.2.1  For  2  <  fc  <  m  ,  uj  =  *‘*^*-^** 

Proof  6k,  the  determinant  of  "Ek,  is  nonzero  and 

Bkul  =  (H-  WfcV*)  u^ 

=  ul  +  u;*v*uj’ 

=  WfcAk-i  -b  WfcVfcuf 

=  w*  (Afc-i  +  v*u*) 

=  uJkAk. 

Lemma  2.2.2  For  2<  k  <m,  AkX  =  b*. 
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Proof  By  definition  of  b*  and  A*,  we  have 

bfc  =  Efcbifc-i 

=  EfcEjt_ibfc_2 
=  EfcEfc_i . . .  E2 Aix 
=  AfcX. 

Lemma  2.2.3  For  2  <  k  <  m,  x*  =  9^^ufkk>k- 

Proof  By  Lemmas  2.2.1  and  2.2.2, 

Xk  =  U^x 

=  ^^‘u>tAjtx 
=  d^^Ukhk- 

Lemma  2.2.4  For  2  <  k  <  m,  b/t_i  =  b/t  —  xjt  v^. 

Proof  Starting  with  Lemma  2.2.2,  we  obtain 

1  ~  Ajt_iX 

=  (A*  -  VkUk)  X 

=  A*x  -  VifcU^x 

=  bfc  -  Vfcx*. 

Lemmas  2.2.3  and  2.2.4  provide  the  same  update  formula  as  Section  1.5  but 
they  show  that  the  components  x*  (m  >  k  >  2)  can  be  computed  along  with 
the  vectors  bt_i  (m  >  k  >  2)  without  additional  work.  By  the  time  bj  is 
computed,  the  only  component  of  x  that  remains  unknown  is  ii.  At  that  stage, 
instead  of  solving  the  whole  triangular  system  AjX  =  bi,  we  only  need  to  solve 
the  first  equation,  of  the  form  AnXi  =  ufbj.  In  summary,  we  obtain  the 
following  method. 

Algorithm  2.2  (to  solve  Ax  =  b  ) 

Let  bm  =  b. 

For  k  =  m  downto  2,  let  Xk  =  9^^Ufk^k  and  bfc_i  =  b/t  —  Xjk  v*. 

Let  xi  =  Af/ufbi- 
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This  algorithm  actually  depends  on  the  auxiliary  vectors  i4fi^  uf  and 
{2  <  k  <  m)  which  are  none  other  than  the  vectors  (1  <  A:  <  m)  (cf. 

Lemma  2.2.1).  The  maximum  number  of  multiplications  required  by  this  method 


E2*  +  1 


+  m. 


k=:m 

This  is  essentially  the  same  as  the  maximum  number  of  multiplications  required 
to  solve  a  system  LUx  =  b  where  L  and  U  are  tri2inguleir  matrices. 

Note  that  only  the  first  k  components  of  bfc  are  needed  for  these  computa¬ 
tions.  Therefore,  for  each  k,  we  may  replace  bt  by  its  projection  b*  onto  the 
subspace  of  generated  by  the  first  k  unit  vectors. 


Example; 

/8  1  6\ 

Consider  the  3x3  system  Ax  =  b  where  A  =  I  3  5  7  I  and  b  = 

\4  9  2/ 

We  know  from  Section  1.4  that  the  factorization  A  =  EjEjAi  entails  the 
following  auxiliary  quantities: 


uf  =  ( 1  ♦  ♦  )  i 

I  *) 

^3  “  (  74  37  2  *  3  1 

A  straightforward  application  of  Algorithm  2.2  yields 


The  final  result  is  the  same  as  in  Section  1.5: 


X  = 


1 

-2 

3 
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2.3  Solution  of  it  A  =  7 

As  in  Section  2.2,  we  can  streamline  the  method  explained  in  Section  1.6  and 
avoid  solving  the  complete  triangular  system  involving  Ai-  For  1  <  ^  <  m,  define 
7^  as  the  projection  of  7  onto  the  subspace  S*  of  generated  by  the  first  fc 

unit  vectors.  Then  define  Tf/t  as  the  unique  solution  of  the  system  =  7^. 

Now,  consider  the  following  Lemmas: 

Lemma  2.3.1  For  \  <  k  <  m  ,  Wk  is  in  St- 

Proof  The  last  m  —  k  equations  of  WjtAfc  =  7;^  constitute  a  full  rank 
triangular  system  with  null  right  hand  side.  Its  solution,  a  vector  made  of  the  last 
m  —  k  components  of  Tf^,  must  be  zero. 

Lemma  2.3.2  For  2  <  A:  <  m,  (Wk  —  ^k-i)  At  =  (7*  —  Ifk-iVk)  , 
where  7^  is  the  A:'*  component  of  7  or,  equivalently,  of  7^  . 

Proof  Using  essentially  the  definitions  of  7^  and  7fc_i,  we  obtain 

WfcAfc  =  7fc 

=  7fc-i  +  Tfcul 

=  Wk-iAk-i  +  TfcUfc 

=  Wk-i  (Ak  -  VkuJ)  +  jk^J 

=  Ifk-iAk  +  (7fc  -  u^. 

Lemma  2.3.3  For  2  <  k  <  m,  TTk  =  Hk-i  +  &k^  (7fc  —  Wfc. 

Proof  The  matrix  Afc  is  non-singular  and  we  have 
(Iffc  -  TCfc-i)  Afc  =  (7fc  -  TTfc-iVfc)  Ufc’ 

=  (7*  —  ^*-iVfc)  9'^^uJkA.k  ■ 

By  Lemma  2.3.1,  only  the  first  component  of  Iti  can  be  nonzero.  Therefore, 
the  system  WiAi  =  7i  can  be  solved  in  one  scalar  division  as  opposed  to  the 
system  ffiAi  =7  which  requires  a  triangular  matrix  division.  Then,  Lemma 

2.3.3  shows  that  the  sequence  Wfc  {2  <  k  <  m)  cem  be  computed  just  as  easily 
as  the  sequence  TTk  {2  <  k  <  m)  of  Section  1.5.  Finally,  since  7  =  7^,  we  have 
JT  =  Wm.  In  summary,  we  obtain  the  following  method. 
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Algorithm  2.3  (to  solve  ttA  =  -y  ) 

_  .IT* 


Let  ffi  =  7i  A^i  u[. 


For  A:  =  2  to  m  ,  let  =  yk  —  T^k-i'Vk  and  iTk  =  ^ 

Let  fji  • 


Once  again,  this  algorithm  actually  depends  on  the  vectors  u^Aj^ ^  (1  <  A;  < 
m).  The  maximum  number  of  multiplications  required  is 

m 

1  +  ^  2A:  =s  +  m. 


Example; 

/8  1  6\ 

Consider  the  3x3  system  ir  A  =  7  where  A  =  I  3  5  7  j  and  7  = 

V4  9  2/ 

(6  0  —6).  We  know  from  Section  1.4  that  the  factorization  A  =  E3E2A1 

entails  the  following  auxiliary  quantities: 


uf  =  (  1  ♦ 

*) 

«|» 

II 

7 ::: 

^2  =  {-^ 

\  ♦) 

Q-l  _  40 
^2  —  37 

U>3  =  (  ^ 

34  1  ^ 

37  2  / 

«  , 

II 

1 

A  straightforward  application  of  Algorithm  2.3  yields 
iri=6i(l  *  ♦)  =  (!  *  *) 


6  =  0  -  (  f  *  ♦ ) 


^2  =  (!  *  *)  -  V^i-ro  i  *)  =  (I?  *) 

«.  =  -6  -  (!?  -if  •)M  = 

»3  =  (I?  -i,  •)  +  ^^  (^  -U  =  (1  -2  1) 

The  final  result  is  the  same  as  in  Section  1.6: 


•  =  -f 


IT  =  (1  -2  1). 
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2.4  Computation  of  the  Factors  E*  (2  <  i  <  m) 

The  system  u>jfcAfc-i  =  u^  defining  tJt  can  be  written 

/A  0  0  \ 

iJk  I  p  oc  o|=(0  1  0) 

\B  c  D/ 

where  A  is  the  {k  —  1)  x  (A:  —  1)  leading  submatrix  of  A  or  At_i,  and  a  is  the 
diagonal  element  of  A  or  Ak-i-  Assuming  that  At_i  is  non-singular,  this  system 
is  equivalent  to 

u)k  =  (  TT  1  0)  and  irA  = -p. 

From  the  relationship  Ak-i  =  (I  +  Vk-i^k-i)  (I  +  Vk-2<^k-2)  • . .  +  V2W2) 

Ai ,  we  derive  A  =  (I  -I- )  (I  +  Vfc_2U>fc-2)  •  ■ .  (I  +  V2u;2)  Ai  where  v/  is 
the  (A:  —  1)  X  1  leading  submatrix  of  vj,  Ui  is  the  1  x  (A:  —  1)  leading  submatrix  of 
u>(  and  Ai  is  the  (A:  —  l)x(A:  —  1)  leading  submatrix  of  Aj. 

When  we  compute  the  vector  we  already  know  the  vectors  a;/  (2  < 
/  <  A:  —  1)  and  hence  the  vectors  (2  <  /  <  A:  —  1).  Therefore,  we  can  apply 
Algorithm  2.3  to  solve  the  (A:  —  1  )-dimensional  system  ttA  =  — p. 

The  mzLximum  number  of  multiplications  or  divisions  required  to  compute 
(Uki^k)  by  this  method  is 

k^  —  k  +  k  -h  k  —  I  ^  k^  +  k 

5md  the  maximum  number  of  multiplications  required  to  factorize  A  is 

m 

fc=2 

This  is  essentially  the  same  as  the  maximum  number  of  multiplications  required 
to  compute  the  triemgulzir  factorization  A  =  LU  by  Gaussian  elimination. 
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Example: 

/8  1  6\ 

Consider  the  3x3  matrix  A  =  [3  5  7  I  . 

\4  9  2/ 

•  We  have  An  =8. 

•  The  factor  E2  =  I  +  V2u;2  is  computed  from  0^2  A 1  =  uj. 

Q  =  5 


rr(8) 

=  ( 

-3 

) 

JT 

=  ( 

3 

8 

) 

UJ2 

_  1 
5 

(- 

3 
'  8 

1 

0  =  i-T. 

i  0 

62 

=  1 

+ 

( 

_3. 

40 

i  *) hi 

_  37 

40 

•  The  factor  E3  =  I  4-  vsWs  is  computed  from  u>3A2  =  uj. 

=  (-1  *) 


i-k  *)  + 

-9  - 

*)[: 

M  lOf -A. 

^  1  37  V  40  5  / 

II 

1 

^*^3  ~  2  (  37 

68 

37 

1)  =  in 

34  1  \ 

37  2  ' 

».  =  1  +  =  -m 


2.5  Simplification  Applicable  to  Non-Spikes 

We  saw  in  Section  1.7  that  if  column  /  is  not  a  spike  of  A,  the  Basic  Method 
does  not  require  the  computation  of  the  auxiliwy  vector  u;j.  This  simplification  is 
not  directly  transferable  in  the  Streamlined  Method,  because  algorithms  2.2  and 
2.3  appear  to  require  the  computation  of  ail  the  intermediate  quantities  like  x/ 
or  Wf,  smd  hence  of  all  the  vectors  W/. 
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However,  if  the  spikes  of  A  were  its  rightmost  columns,  the  leading  t  x  t 
submatrix  T  sissociated  with  the  t  non-spike  columns  of  A  would  be  trizuiguiar. 
Then,  the  vectors  (I  <  I  <  t)  would  represent  the  rows  of  T“*.  Moreover, 

the  steps  of  algorithms  2.2  and  2.3  corresponding  to  spike  columns  k  could  be 
carried  out  as  earlier  with  the  vectors  utjt,  while  those  corresponding  to  nonspike 
columns  I  could  be  replaced  by  solving  a  system  explicitly  involving  the  triangular 
matrix  T  instead  of  the  vectors  u»/. 

In  order  to  implement  this  idea  when  the  matrix  A  is  arbitrary,  we  can  per¬ 
mute  the  columns  of  A  according  to  the  permutation  matrix  Q  that  sends  the 
spikes  of  A  to  the  right  while  preserving  the  order  of  the  non-spikes  and  that  of 
the  spikes,  and  permute  the  rows  of  A  symmetrically,  i.e.  according  to  the  per¬ 
mutation  matrix  Q^.  Let  A'  =  Q^AQ  be  the  resulting  matrix.  Let  q  be  the 
permutation  of  {l,2,...,m}  induced  by  Q  (i.e.  Ut^k)  =  Q^Ufc)-  Then  we  have 
the  following  Lemma: 

Lemma  2.5.1  Column  j  is  a  spike  of  A  if  and  only  if  column  q(j’)  is  a  spike 
of  A'. 

Proof  If  column  j  is  a  spike,  then  3i  i  <  j  and  Aij  ^  0  .  By  construction, 
if  column  j  is  a  spike  of  A,  then  i  <  j  q(i)  <  q(j)  .  In  addition, 

Q^AQ  =  ufAuj  =  A,j  ^  0  .  Therefore, 

column  q(y )  is  a  spike  of  A'. 

Conversely,  if  column  q(j)  is  a  spike  of  A',  then  3i  q(i)  <  q(j)andAq( ^  0 
.  Either  i  <  j  or  j  <  i.  In  the  former  case,  3i  i  <  j  and  ^  0  .  In  the  latter 
case,  column  j  which  stood  to  the  left  of  column  i  in  A  will  steuid  to  the  g 
column  :  in  A'.  Either  way,  column  j  must  be  a  spike  of  A. 


By  Lemma  2.5.1,  A  and  A'  have  the  same  number  s  of  spikes  and  the  same 
number  t  =  m  —  s  of  non-spikes.  Since  the  i  leftmost  columns  of  A'  are  not  spikes, 
the  decomposition 


A'  =  A'l  -I-  VjuJ  +  VjuJ’ 


+  •••  +  v'  u 


_  T^/ 

m**Tn  ‘^m  ^m  — 1 


E'  a; 
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I.  In 


satisfies  Vj  =  Vj  =  ■  •  ■  =  vj  =  0  and  E2  =  Ej  =  •  •  ■  =  E',  = 
particular,  AJ  is  triangular. 

To  solve  Ax  =  b,  we  consider  the  equivalent  system  A'x'  =  b'  where 
x'  =  Q^x  and  b'  =  Q^b.  The  last  s  components  of  x'  are  computed  as  in 
Section  2.2  using  the  vectors  (m  >  k  >  t),  and  the  first  t  components  of  x' 
are  solved  through  the  first  t  equations  of  the  triangular  system  AJx'  =  bj. 

To  solve  IT  A  =  'Y  ,  we  consider  the  equivalent  system  tt'A'  =  ■y'  where 
tt'  =  jtQ  cmd  y'  =  yQ  .  The  vector  IrJ  is  given  by  the  triangular  system 
^',A(  =  y't  and  the  last  s  vectors  W'/^  are  computed  as  in  Section  2.3  using  the 
vectors  {i  <  k  <  m). 

To  compute  the  non-trivial  factors  E'^  (t  <  k  <  m),  we  decompose  the  system 
<*;^A^_j  =  into  =  a~^  ( ff  1  0)  and  irA  =  —p  as  shown  in 

Section  2.4  and  solve  the  latter  system  as  indicated  in  the  previous  paragraph. 

That  way,  all  computations  can  be  carried  out  solely  with  the  auxiliary  quan¬ 
tities  (t  <  k  <  m)  corresponding  to  spikes  (v'^  ^  0).  Because  the  spikes 

have  been  shifted  to  the  right  in  A',  the  number  of  nonzeros  above  the  diagonal 
in  the  spikes  may  increase.  However,  the  number  of  nonzeros  in  the  vectors 
remains  unchanged  as  the  following  Lemma  indicates: 

Lemma  2.5.2  If  column  k  is  a  spike  of  A,  then  and 

■ 

Proof  Let  k'  =  q(fc).  We  have,  by  definition  of  , 

and,  by  definition  of  utk, 

w*A*_i  =  u* 

Q^Ak-iQ  =  ulq  =  ujfc)  =  u[, 

Partition  Q^Afc_iQ  into  (y  t)  "K  is  a  k^  X  k'  submatrix.  Then 

Z  =  0,  X  and  T  have  full  rank  and  X  equals  the  leading  k'  x  k'  submatrix  of 


A'  =  Q^AQ.  Therefore,  the  last  m  —  k'  components  of  0;^,  and  of  are 

zeros  while  their  first  k'  components  are  solutions  of  the  same  well  determined 
system.  Thus  ~  ■ 

Finally,  since  the  first  k'  components  of  and  Qj'Vk  are  equal,  we  have 

^q(fc)  =  1  +  ‘*^q(ifc)Vq(fc) 

=  1  +  wjtQfcQfcVifc 

=  1  +  UJfcVfc 

=  Ok . 


Example: 

This  example  shows  how  to  permute  and  factorize  a  5  x  5  matrix  whose 
columns  2  zmd  4  are  not  spikes,  and  how  to  solve  systems  like  Ax  =  b  or  ttA  =  y. 
In  this  case,  the  permutation  Q  will  simply  interchange  columns  3  and  4. 


/ 1 

0 

1 

0 

(^\ 

0 

1 

0 

0 

0  ' 

1 

-1 

0 

1 

0 

1 

b  = 

6 

0 

1 

1 

1 

6 

\  1 

1 

0 

1 

1) 

Vs/ 

7  =  (1  3  3  2  3) 


/ 


A'  = 


1 

0 

0 

1 

0 

1 

0 

0 

0 

1 

0 

1 

1 

1 

0 

b'  = 

6 

-1 

0 

0 

1 

1 

6 

1 

1 

1 

0 

1/ 

Vs/ 

V  =  (l  3  2  3  3) 
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a;  =  a'2  =  A'  = 


/ 1 

* 

♦ 

* 

0 

1 

* 

* 

♦ 

0 

11 

> 

=► 

E2  =  I 

0 

1 

1 

* 

♦ 

-1 

0 

♦ 

1 

♦ 

0 

11 

> 

E'  =I 

V  1 

1 

1 

0 

1/ 

/ 1 

* 

* 

1 

0 

0 

1 

* 

0 

♦ 

V4  = 

1 

0 

1 

1 

1 

♦ 

* 

-1 

0 

♦ 

1 

* 

\  1 

1 

1 

0 

1) 

/ 1 

♦ 

* 

1 

i  1 

0 

0 

1 

♦ 

0 

0 

< 

11 

0 

0 

1 

1 

1 

0 

1 

-1 

0 

♦ 

1 

1 

V  1 

1 

1 

0 

1> 

' 

•  The  factor  E4  =  I  +  is  computed  from  w^Aj  =  uj . 


^(J  ;)  =  (1  0) 

TT  =  (1  0) 


1-^(1  0  *  1 
1  +  (1  0  * 


*) 


*) 


=  (1 

/i\ 


0 

1 


0  ♦ 


=  2 


♦) 


•  The  factor  E5  =  I  4-  VgUij  is  computed  from  utjAJ,  =  uj. 


IT 


♦  * 
1  * 
1  1 
0  ♦ 


=  1 

=  (-1  -1  -1  0) 


W3  =  i-l  0  -1  *) 
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•  The  vector  x  is  computed  from  E5E4A3X'  =  b'. 


•  The  vector  tt  is  computed  from  ir'E5E4  A3  =  “i' . 

W'A;  =(132**) 

^3  =  (1  1  2  *  *) 
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«0  «D  00 


w;  =(i  1  2 


*  * )  + 


CHAPTER  3.  NUMERICAL  STABILITY 


3.1  Introduction 

In  Chapter  2,  we  showed  how  to  factorize  a  matrix  A  under  the  assumption 
that  the  matrices  At  (1  <  A:  <  m)  were  non-singular.  In  this  Chapter,  we 
decompose  this  assumption  into  two  conditions  and  show  that  a  rescaled  version 
of  the  same  factorization  can  be  obtained  under  only  the  first  condition.  We 
recall  that  if  A  is  non-singular  then  its  columns  can  be  permuted  so  that  the 
resulting  matrix  satisfies  the  first  condition,  and  indicate  a  procedure  to  carry 
out  simultaneously  the  permutation  and  the  factorization.  This  procedure  offers 
a  numerical  stability  similar  to  that  of  LU  decomposition  with  partial  pivoting. 


3.2  Definitions  and  Notation 

A  unit  upper  triangular  matrix  is  an  upper  trianguleu’  matrix  whose  entries 
on  the  diagonal  axe  equal  to  one. 

Let  A  be  a  squaure  matrix.  Let  B  be  a  square  submatrbt  of  A.  We  define  the 
minor  of  A  associated  with  B  as  the  determinant  of  B.  A  leading  minor  of  A  is 
a  minor  associated  with  a  leading  submatrix  of  A. 

Let  A,v  and  (jJ  be  respectively  an  m  x  n,  an  m  x  1  and  a  1  x  n  matrix.  Then, 
for  ib  <  m  and  k  <  n,  (A)fc,  (v)*  and  (u;)fc  denote  respectively  the  k  x  k  leading 
submatrix  of  A,  the  fc  x  1  leading  submatrix  of  v  and  the  1  x  fc  leading  submatrix 
of  uf.  Unless  indicated  otherwise,  A*,  v*  md  u;*  are  abbreviations  for  (Afc)fc, 
(Vfc)t  and  {u}k)k- 


3.3  Conditions  on  the  Matrices  Afc  (1  <  A:  <  m) 

In  Chapter  2,  we  assumed  that  the  matrices  Afc  (1  <  fc  <  m)  were  non- 
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singul^l^.  Because  of  their  structure,  we  have 

m  m 

det  A*  =  det  Afc  Au  =  det  At  n 

l=k+l  l=k+l 

Therefore,  the  matrices  A/t  (1  <  ^  <  m)  are  non-singular  if  and  only  if  the 
matrices  A*  (1  <  <  m)  and  the  matrix  Ai  are  non-singular. 

It  is  well  known  that  if  A  is  non-singular,  then  its  columns  can  be  permuted 
according  to  a  permutation  matrix  Q  such  that  the  leading  square  submatrices  of 
AQ  are  non-singular.  Actually,  the  columns  of  A  can  be  selected  step  by  step:  at 
step  k,  a  column  is  chosen  for  position  so  as  to  maximize  the  absolute  value  of 
the  resulting  A:**  leading  minor.  In  practice,  this  procedure  yields  well  conditioned 
leading  submatrices,  as  observed  in  the  LU  decomposition  zJgorithm  with  column 
interchanges.  Therefore,  we  shall  follow  a  similar  procedure. 

Regarding  the  non-singularity  of  (AQ)i,  i.e.  the  absence  of  zeros  on  the 
diagonal  of  AQ,  we  shall  see  in  Section  3.6  that  the  issue  is  rendered  moot  by 
rescaling  the  vectors  representing  the  factorization. 


3.4  LU  Decomposition  with  Column  Interchanges 

In  this  Section,  we  recall  some  useful  properties  of  the  LU  decomposition  of 
square  matrices  (Murty,  1976). 

Lemma  3.4.1  For  any  m  x  m  matrix  A,  there  exist  a  lower  triangular 
matrix  L,  m  —  1  transposition  matrices  T*  (1  <  A:  <  m)  and  m  —  1  unit  upper 
tri2mgular  matrices  Ut(l<A:<m)  such  that 

L  =  AT,UxT2U2...T^_iU„,_,. 

Proof  (and  algorithm  3.4)  Let  A^®^  =  A.  Given  let  jk  satisfy 

If  =  0,  is  singular.  Let  Tt  =  Ufc  =  I. 

Otherwise,  interchange  columns  k  and  jk,  i-e.  postmultiply  by  the 

transposition  matrix  T/tj*  (abbreviated  to  Tt)  derived  from  the  identity  matrix 
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by  interchanging  columns  k  and  jk.  Then  zero  out  the  entries  of  row  k  to  the  right 
of  column  k  by  adding  appropriate  multiples  of  column  k  to  columns  +  1, . . . ,  m, 
i.e.  postmultiply  by  a  unit  upper  triangular  matrix  U*  with  the 

appropriate  entries  in  row  k. 

In  either  case,  we  have 

Xii^)  ^  A^^-'^^TkUk  for  1  <  Jb  <  m  -  1 

where  the  strictly  upper  triangular  part  of  has  zero  entries  in  its  k  first  rows. 
After  m  —  1  iterations,  we  obtain  =  A^®^TiUiT2U2  . . . 

where  is  lower  triangular. 

Lemma  3.4.2  For  0  <  A:  <  m  -  1,  let  A^*^)  =  AT1U1T2U2  . . .  T^Ut 
be  the  A:'*  iterate  in  the  LU  decomposition  of  A.  Then  A^*'^  =  AT1T2  . . .  TtU'^ 
where  U'^  is  a  unit  upper  triangular  matrix  whose  last  m  —  k  rows  equal  those 
of  the  identity  matrix. 

Proof  (by  induction)  A^°^  =  AI.  Assume  that 
A(fc-i)  ^  ATiT2...Tfc_iU;_i  where 
Then 

;)(!  ?) 

/V  QS\ 

-  Vo  s ; 

fl  0\  /V  QS\ 

Vo  s; Vo  1  ) 

= 

and 

A(*)  =  A^*-^>TkUfc 

=  AT,T2...Tfc_iUUTfcUfc 
=  ATiT2...T*_iT*V;_iUfc 
=  ATiT2...T*_iTfcU; 
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where  V^_j,  Ufc  and  UJt  =  are  unit  upper  triangular  matrices  whose  last 

m  —  k  rows  equal  those  of  the  identity  matrix. 


Corollary  3.4.3  There  exists  a  lower  tri2ingular  matrix  L,  a  unit  upper  tri¬ 
angular  matrix  U  and  a  permutation  matrix  Q  such  that  LU  =  AQ. 


Corollary  3.4.4  For  \  <  I  <  k,  the  leading  minor  of  and  that  of 
AT1T2  . . .  Tfc  are  equal: 


det(A(*))/  =  det(ATiT2...Tfe)/. 


Corollary  3.4.5  If  A  is  non-singular,  the  first  k  leading  minors  of  the 
product  ATi  T2  . . .  Tt  are  nonzero: 


det(ATiT2...Tfc),  =  det(A(*'))j 


—  "^1,1  -^2,2  •••  ^1,1 


= 


^2,2 


...  A 


(m-l) 


Corollary  3.4.6  If  A  is  non-singular,  the  pivot  is  the  ratio  of  the 
leading  minor  over  the  (fc  —  1)**  leading  minor  of  AT1T2  . . .  T^: 


,(fc)  det(A(*))fc  det(ATiT2...Tfc)fc 

^k,k  — 


det(A(*^)fc_i  det(ATiT2  . . .  Tifc)jfc_i 

Therefore,  in  terms  of  the  matrix  AT1T2  •  •  •  Tfc,  the  following  column  selec¬ 
tion  rules  are  equivalent: 

(1)  given  the  first  jfc  —  1  columns,  select  the  column  so  as  to  maximize  the 
absolute  value  of  the  resulting  k*^  pivot. 

(2)  given  the  first  jfc  —  1  columns,  select  the  column  so  as  to  maximize  the 
absolute  value  of  the  resulting  leading  minor. 
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3.5  An  Alternate  Computation  of  the  Pivots 

In  this  Section,  we  show  how  to  compute  the  potential  pivots  through  some 
auxiliary  vectors  <Tk  {I  <  k  <  m),  without  explicitly  carrying  out  any  LU  decom¬ 
position. 


Lemma  3.5.1  Let  A  be  a  non-singular  m  x  m  matrix.  Let  (1  <  ^  <  m) 
be  the  sequence  of  pivots  generated  by  the  LU  decomposition  of  A  as  described 
in  Section  3.4.  Let  Q  =  TiT2...Tm-i  be  the  permutation  matrix  resulting 


from  the  same  LU  decomposition.  Let  A'  =  AQ.  Let 


V  Pk 


-1  afe 

Cuk 


,  where  ot 


is  a  scalar,  represent  the  k  x  k  leading  submatrix  of  AT1T2  . . .  T*,  and  hence  the 
k  X  k  leading  submatrix  of  A'.  Let  iv  =  {1 . . .  k). 

Then,  for  1  <  1;  <  m,  the  m-vector  (Tk  =  (  — pt  A^t_i  1  0 )  is  well 

defined  and  the  pivot  l3k  is  given  by  the  scalar  product  /3k  =  "^k  ^Kk- 

Proof  By  Corollary  3.4.5,  A'k-i  is  non-singular  and  by  Corollary  3.4.6, 
we  have 

det 


I3k  = 


A'fc_i  afc 
Pk  Oik 


det  A' k-i 
-1 


=  ttfc  —  PfcA'jt-iafc 

=  <Tk  A'f^k- 

Note  that  for  l<A:<m,  because  A  is  non-singular  (cf.  Corollary 

3.4.3). 


The  algebraic  relationships  mentioned  in  Section  3.4  hold  independently  of 
the  pivot  selection  rule.  Therefore,  at  step  k,  all  potential  pivots  can  be  computed 
as  products  of  the  form  where  c  is  one  of  the  unselected  columns  of  A/c,,  i.e. 
one  of  the  m  —  fc  -t-  1  rightmost  columns  of  A'^,.  In  other  words,  the  potential 
pivots  are  the  m  —  fc  -f  1  rightmost  elements  of  the  vector  A  =  WkA'  —  where 

t\  i\ 

K  = 
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3.6  Factorization  with  Column  Interchanges 

If  A  is  non-singvilar,  the  column  interchange  method  described  above  provides 
a  permutation  matrix  Q  =  T1T2  . .  .Tm  such  that  A'  =  AQ  has  nonzero  leading 
principal  minors.  The  following  theorem  explains  why  the  factorization  and  the 
column  permutation  cam  be  carried  out  simultaneously. 

Theorem  3.6.1  Let  A  be  a  non-singular  m  x  m  matrix.  Let  Q,  A'  and 
(<rk,0k)  (l<fc  <  m)  be  as  defined  in  Lemma  3.5.1.  Let  ajt  (1  <  fc  <  m)  denote 
the  diagonal  elements  of  A'.  If,  for  1  <  A:  <  m,  a*  ^  0,  then  A'  can  be  factorized 
ais  described  in  Chapter  2.  and  the  sequence  (ufk,^k)  (2  <  A:  <  m)  resulting  from 
this  factorization  satisfies 


u>fc  =  and  6k  =  oik^l^k  for  2  <  A:  <  m  . 

Proof  The  matrix  A*.,  is  non-singular  and  we  have 


<rkK-i  =  (-pkA'tli  1  0) 

Therefore  <r*Aj_j  = 


6k  =  det  E* 


det  A'k 
det  A^_j 


det 


det 


A'k-i  0  0\ 

Pk  0  I  =  (0  Ok  0)  = 

♦  *  *  / 


and  w*  =  .  Finally,  Wc  have 

/l^k-i  a*  N  _ ^ 

\  Pk  o^k )  _  at  -  PkA'k-i^k 

fUk-i  0  \  “  at 

\  Pk  oik) 


T 

at  ui. 


h 

at' 


Since  <Tt  can  be  regarded  as  the  value  that  ufk  would  teJce  if  at  were  equal  to 
one,  it  can  be  computed  at  least  as  easily  €is  Wf  As  a  matter  of  fact,  the  method 
used  in  Section  2.4  to  compute  ufk  entails  the  computation  of  <Tt  =  (  1  0 ) 

where  ?rAt-i  =  — pt  is  solved  by  Algorithm  2.3. 

In  addition,  Theorem  3.6.1  implies  that,  for  2  <  k  <  m,  6'j^^uJk  =  ^^^(Tk- 
Therefore,  in  Algorithms  2.2  and  2.3,  the  sequence  {(Tki^k)  (2  <  A:  <  m)  can 
replace  the  sequence  (ufk,6k)  (2  <  k  <  m)  which  need  not  be  computed. 
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Finally,  if  A  is  non-singular  then  the  sequence  {(Tk^^k)  (1  <  fc  <  m)  is  well 
defined  whereas  the  existence  of  (2  <  t  <  m)  also  depends  on  whether 

the  coefficients  are  different  from  zero. 

All  these  considerations  lead  us  to  redefine  the  factorization  of  A  in  terms  of 
the  sequence  ((Tk.i^k)  {I  <  <  m)-  Using  notation  of  Section  3.5,  we  end  up 

with  the  following  method. 

Algorithm  3.6  (to  generate  Q  and  factorize  A'  =  AQ) 

Let  A'  =  A. 

For  k  =  1  to  m,  do  the  following: 

•  Solve  jrA'jt-i  =  ~Pk- 

•  Let  tTfc  =  (  JT  1  0 ). 

•  Compute  the  vector  of  potential  pivots  A  = 

•  Select  a  pivot  column  such  that  lAj^j  =  maxfc<)<m  lAji|. 

•  Let  A'  =  A'Tkjk- 

•  Let  l3k  = 

•  U  0k  =  0  then  stop  (A  is  singular). 


Example:  Consider  the  3  x  3  matrix  introduced  in  Section  1.2: 


/8  1  6' 

A  =  3  5  7  1  . 

\4  9  2^ 

•  Step  1  Column  1  has  the  entry  of  largest  absolute  value  in  row  1.  There¬ 
fore,  column  1  remadns  in  first  position. 

<Ti  =  (1  ♦  *) 

A  =  (1)  (8  1  6)  =  (8  1  6) 

01  =  S 

•  Step  2  We  compute  the  vector  <72: 


<T2 


a  ;)  = 


(0  1) 
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flTj  =  (  «■  1  *  )  where  ir  ( 8  )  =  (  -3  ) 

^2  =  (-1  1  *) 


The  potential  pivots  are  given  by 


A  =  (-t  1)  0  ?)  =  T)- 

We  select  02  -  T  -  ^  column  3  which  moves  into  second  position. 
•  Step  3  We  compute  the  vector  (Ti'. 

/8  6  0\ 

(Til  3  7  0]=(0  0  1) 

\4  2  ij 


fl-3  z=  ( 1 )  where 


®)  =(-4  -2) 


=  (-5  *) 


^2 


=  i-h  *)  + 


-  (4  •)(!)] 


^3  —  (~T9  19 


The  only  remaining  pivot  is  given  by 


A  =  (-il  ^  1) 


180 

19 


a  -  180 

03  —  19 


34 


<T,  =  (  1 


♦  *)  ^1=8 


8  *  * 

3  7* 

4  2  9 

8  6* 

3  7* 

4  2  9 

5  6  1 

3  7  5 

4  2  9 


^2=(-|  1  *)  /?2  =  -^ 
^1)  =  W 


3.7  Rescaling  of  the  Matrix  A 

To  reduce  the  number  of  divisions  during  the  factorization  of  a  matrix,  it  may 
be  worthwhile  to  rescale  the  columns  of  the  matrix  so  that  the  resulting  pivots 
(Jfc  (1  <  <  ^2)  are  equal  to  one.  For  this  resc2iling  to  speed  up  the  factorization, 

the  number  nz  of  nonzero  elements  off  the  diagonal  must  be  less  than  the  number 
of  divisions  by  I3k  : 

^  m(m  -  1) 

nz  <  ^(m  ~k)  = - - - . 

k=l 

Therefore,  the  density  of  the  matrix  must  be  less  than  | .  This  procedure  does  not 
affect  the  number  of  divisions  involved  in  solving  the  system  Ax  =  b  for  instance, 
because  of  the  unsealing  of  x. 

Example: 

/8  1  6\ 

Consider  the  3x3  matrix  introduced  in  Section  1.2  A  =  I  3  5  7  I  . 

V4  9  2/ 

After  the  example  of  Section  3.6,  it  is  easy  to  see  that  the  matrix 


h2is  all  its  pivots  equal  to  one. 


<Ti  =  (  1  ♦  *  ) 

^2  =  (-|  1  *) 
^3  =  (  — H  -jg  1  ) 
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CHAPTER  4.  SPIKE  REDUCTION 


4.1  Introduction 

The  method  presented  in  Section  2.5  to  factorize  a  matrix  A  and  to  solve  a 
system  Ax  =  b  or  ir  A  =  'y  will  work  better  if  the  number  of  auxilizirj’  vectors  cjfc 
(or  (Tjt)  is  smaller,  and,  given  their  number,  if  their  size  is  smaller. 

The  most  obvious  way  to  achieve  these  goals  is  to  reduce  the  number  of  spikes 
of  the  matrix  A  and,  whenever  possible,  to  shift  those  spikes  towards  the  left.  To 
that  end,  we  describe  three  myopic  algorithms,  inspired  by  Hellerman  and  Rarick’s 
procedure  (Hellerman  and  Rarick,  1972),  that  reorder  the  matrix  by  selecting 
some  rows  and  columns,  assigning  them  a  position,  deleting  them  and  repeating 
the  process  on  the  resulting  submatrix  imtil  all  rows  or  all  columns  have  been 
assigned. 


4.2  Definitions  and  Notation 

Let  i  and  j  denote  a  row  and  a  column  index  of  A.  If  A ^  0,  we  say  that  row 
i  “intersects”  column  j,  that  column  j  “intersects”  row  i  or  that  row  i  eind  column 
j  “intersect” .  We  define  the  following: 


ROWSPAN(i) 

COLSPAN(i) 

ROWCOUNT(i) 

COLCOUNTO) 

ROWTALLY(i) 

COLTALLYO) 


the  set  of  columns  j  intersecting  row  i. 

the  set  of  rows  i  intersecting  column  j . 

the  number  of  nonzero  entries  in  row  i. 

the  number  of  nonzero  entries  in  column  j. 

the  ninnber  of  nonzero  entries  in  the  columns  intersecting  row  i. 

the  number  of  nonzero  entries  in  the  rows  intersecting  column  j. 
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ROWSPAN(j)  =  {j  :  Aij  ^  0} 
COLSPAN(j)  =  {i  Aij  ^  0} 


ROWCOUNT(i)  =  1rOWSPAN(01 
COLCOUNT(j)  =  |COLSPANO)| 


A 

ROWTALLY(0  = 


COLCOUNTO) 

j-.Aij  #0 


A 

COLTALLY(i)  = 


ROWCOUNT(i) 


A  line  of  a  matrix  A  is  a  row  or  a  column  of  A.  At  iteration  fc,  the  “active 
submatrix”  is  the  submatrix  obtained  after  deletion  of  the  lines  selected  during 
iterations  1, . . . ,  t  —  1. 


4.3  Top>Left  Spike  Reduction  Algorithm 

This  algorithm  essentially  keeps  selecting  from  the  active  submatrix  a  row  to 
be  placed  at  the  top  of  the  unassigned  rows  and  a  matching  column  to  be  placed 
to  the  left  of  the  unassigned  columns.  This  amounts  to  identifying  the  coefficient 
of  the  active  submatrix  to  be  placed  in  the  top-left  comer. 

More  precisely,  at  iteration  fc,  a  row  with  the  fewest  nonzero  entries  is  se¬ 
lected  from  the  active  submatrix  and  assigned  position  fc  (the  highest  available). 
If  some  columns  of  the  active  submatrix  intersect  that  row,  one  of  them  is  as¬ 
signed  position  fc  (the  leftmost  available)  and  the  others  are  sent  to  a  spike-index 
queue.  Otherwise,  a  coliunn  is  removed  from  the  spike-index  queue  according  to 
the  First-In-First-Out  (FIFO)  priority  mle  and  assigned  the  position  fc.  Finally, 
the  active  submatrix  is  updated  by  deletion  of  the  selected  row  and  of  the  columns 
intersecting  it. 

Termination  occurs  when  all  columns  have  been  deleted.  Then  the  undeleted 
rows  are  assigned  the  bottom  positions,  the  columns  remaining  in  the  spike-index 
queue  are  removed  in  FIFO  order  and  assigned  the  rightmost  positions.  Under 
this  algorithm,  deletion  of  all  rows  cannot  occur  before  deletion  of  all  columns. 
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In  addition,  note  that  the  spikes  added  to  the  queue  during  iteration  k  will 
have  zero  entries  above  row  k  and  a  nonzero  entry  in  row  k  of  the  current  matrix, 
and  that  the  only  unassigned  columns  with  nonzero  elements  above  row  k  axe 
already  in  the  spike  queue.  Therefore,  if  no  column  intersects  the  selected  row 
and  no  spike  is  available  in  the  queue  at  iteration  k,  then  the  first  k  rows  of 
the  reordered  matrix  contain  only  fc  —  1  nonzero  columns  and  the  matrix  A  is 
structurally  singular  (i.e.  singular  for  any  values  given  to  the  nonzero  coefficients). 

Finally,  the  FIFO  priority  rule  used  in  removing  the  indices  from  the  spike- 
index  queue  produces  spikes  in  order  of  non-increasing  heights,  which  prevents 
the  structural  singularity  of  the  leading  submatrices,  unless  the  whole  matrix  A 
is  itself  structurally  singular. 

The  Top- Left  spike  reduction  algorithm  is  listed  below: 

Let  A:  =  0. 

Repeat 

Let  k  =  k  + 

T0PLEFT(A:); 

until  all  columns  are  deleted. 

Assign  undeleted  rows  positions  { A:  -H  1 . . .  m } . 

Remove  columns  remaining  in  spike-index  queue  (FIFO). 

Assign  these  columns  positions  {A:  -f-  1 . . .  m}. 

The  procedvire  topleft(A:)  consists  of  the  following  instructions: 

•  Select  from  the  active  submatrix  a  row  minimizing  ROWCOUNT(t).  Brealc 
ties  by  maximizing  ROWTALLY(i)- 

•  Assign  row  i*  position  k. 

•  If  R0WC0UNT(i*)  ^  0,  then 

select  a  pivot  column  jk  from  the  columns  intersecting  row  it; 

eissign  column  jk  position  A:; 

add  the  other  columns  intersecting  row  to  the  spike- index  queue. 
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•  If  ROWCOUNT(ifc)  =  0,  then 

remove  the  FIFO  colunm  from  the  spike  queue  and  assign  it  position  k. 

•  Delete  row  i*  and  the  columns  intersecting  row  I'fc. 


4.4  Bottom-Right  Spike  Reduction  Algorithm 

Instead  of  building  the  new  matrix  from  the  top-left  corner  to  the  bottom- 
right  one,  the  second  algorithm  does  it  in  reverse  direction,  skipping  some  spaces 
for  the  spikes.  It  essentially  keeps  selecting  from  the  active  submatrix  a  column  to 
be  placed  to  the  right  of  the  unassigned  columns  and  a  matching  row  to  be  placed 
at  the  bottom  of  the  unassigned  rows.  This  amotmts  to  identifying  the  coefficient 
of  the  active  submatrix  to  be  placed  in  the  bottom-right  corner. 

At  each  iteration,  a  colunm  with  the  fewest  nonzero  entries  is  selected  from 
the  active  submatrix.  If  some  rows  of  the  submatrix  intersect  that  colunm,  they 
are  assigned  the  positions  at  the  bottom  of  the  submatrix,  the  selected  column 
is  assigned  the  rightmost  position  p  allowing  it  to  fit  on  and  below  the  diagonal, 
and  the  imeissigned  positions  to  the  right  of  p  are  sent  to  a  spike-position  queue. 
Otherwise,  a  position  is  removed  from  the  spike-position  queue  according  to  the 
First-In-First-Out  (FIFO)  priority  rule  and  the  selected  column  is  assigned  that 
position.  Finally,  the  active  submatrix  is  updated  by  deletion  of  the  selected 
column  and  of  the  rows  intersecting  it. 

Termination  occurs  when  all  rows  have  been  deleted.  Then  the  undeleted 
columns  are  assigned  the  positions  remaining  in  the  spike-position  queue.  Note 
that,  under  this  algorithm,  deletion  of  all  coliunns  cannot  occur  before  deletion  of 
all  rows. 

If  no  row  intersects  the  selected  column  and  no  spike-position  is  available  in 
the  queue  at  iteration  k,  then  some  I  columns  of  A  contain  at  most  1—1  nonzero 
rows  and  the  matrix  A  is  structurally  singular. 

The  Bottom-Right  algorithm  obtains  the  spikes  in  order  of  non-decreasing 
heights,  but  it  allocates  them  from  right  to  left.  However,  the  Bottom-Right 
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algorithm,  listed  below,  usually  yields  fewer  spikes  than  the  Top-Left  algorithm; 
Let  k  =  m  +  1. 

Let  pk  =07  +  1. 

Repeat 

Let  k  =  fc  —  1; 
bottomright(A:); 
until  all  rows  are  deleted. 

Assign  the  undeleted  colunms  the  spike-positions  remaining  in  queue. 
The  procedure  bottomright(A;)  consists  of  the  following  instructions: 

•  Select  from  the  active  submatrix  a  column  jk  minimizing  coLCOUNT(j);  break 
eventual  tie  by  maximizing  coltally(j). 

•  Update  the  rightmost  non-spike  position  pk  =  Pk+i  —  colcount(>). 

•  If  coLcouNT(jt)  ^  0,  then 

assign  column  jk  position  pk] 

select  a  pivot  row  from  the  rows  intersecting  column  jk', 
assign  the  pivot  row  position  pk', 

assign  the  other  rows  intersecting  column  jk  the  positions  in  the  range 
Pk  +  l.  “Pk+i  -  1; 

add  Pit  -I- 1, . . .  ,Pfe-i-i  —  1  to  the  spike-position  queue. 

•  If  coLCOUNT(ik)  =  0,  then 

remove  the  FIFO  position  q  from  the  spike-position  queue; 
assign  column  jk  position  q. 

•  Delete  column  jk  and  the  rows  intersecting  column  j*. 

4.5  Composite  Spike  Reduction  Algorithm 

The  third  algorithm  is  a  combination  of  the  first  two.  At  odd  iterations,  it 
selects  and  deletes  a  column  zuid  the  rows  intersecting  it  (using  the  procedure 
bottomright).  At  even  iterations,  it  selects  and  deletes  a  row  and  the  columns 
intersecting  it  (using  the  procedure  topleft). 
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Termination  occurs  when  all  rows  have  been  deleted.  Then,  the  undeleted 
columns  are  added  to  the  spike-index  queue.  Finally  the  columns  remaining  in  the 
spike-index  queue  are  matched  with  the  positions  remaining  in  the  spike-position 
queue. 

Empirically,  this  algorithm  appears  to  combine  the  speed  of  the  first  one  (i.e. 
it  requires  a  small  number  of  iterations)  and  the  efficiency  of  the  second  one  (i.e. 
it  yields  a  small  number  of  spikes). 

The  composite  spike  reduction  algorithm  is  Usted  below: 

Let  k  =  m  +  1. 

Let  pfc  =  m  -t- 1. 

Let  k'  =  0. 

Repeat 

If  not  ail  columns  are  deleted  then 
Let  fc  =  fc  —  1; 

BOTTOMRIGHT(fc), 

If  not  all  rows  are  deleted  then 
Let  k'  =  k'  +  1; 

TOPLEFT(fc'). 

until  all  rows  are  deleted. 

Add  undeleted  columns  to  the  spike-index  queue. 

Match  remaining  spike-indices  eind  remaining  spike-positions. 


4.6  Examples 

In  this  paragraph,  we  apply  the  three  algorithms  described  above  to  an  8  x  8 
matrix.  Although  the  numerical  vsJues  of  the  nonzero  coefficients  are  not  needed, 
they  Me  represented  for  the  sake  of  consistency  with  the  example  used  in  Chapter 
4.  SI  queue  and  SP  queue  denote  the  spike-index  queue  and  the  spike-position 
queue  respectively.  The  columns  belonging  to  the  SI  queue  are  printed  in  italic. 
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4.6.1  First  algorithm: 


2  3  4  5  6  7  8  count  tally 


count  25353334 


Matrix  0 

Iteration  1 

Row  3  is  selected  and  assigned  position  1. 

Columns  4  and  5  are  selected. 

Column  4  is  assigned  position  1  and  column  5  is  sent  to  the  SI  queue. 
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Iteration  2 

Row  5  is  selected  2md  assigned  position  2. 
Columns  2  and  3  are  selected. 
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Column  3  is  assigned  position  2  and  column  2  joins  column  5  in 


the  SI  queue. 
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Iteration  3 

Row  8  is  selected  and  assigned  position  3. 
Column  8  is  selected  and  assigned  position  3. 
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Matrix  3a 


Iteration  4 

Row  1  is  selected  and  assigned  position  4. 
Column  1  is  selected  and  assigned  position  4. 
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Matrix  4a 


Iteration  5 

Row  4  is  selected  and  assigned  position  5. 

Column  5  is  removed  from  the  bottom  of  the  SI  queue  and  assigned  position  5. 
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Iteration  6 

Row  2  is  selected  and  assigned  position  6. 

Columns  6  and  7  are  selected. 

Column  6  is  assigned  position  6  smd  colunm  7  joins  column  2  in  the  SI  queue. 
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Matrix  6a 

After  6  iterations,  all  columns  have  been  deleted.  The  remaining  rows  6  and  7 
are  assigned  the  remaining  positions  7  and  8,  as  are  the  columns  2  and  7  left  in  the 
SI  queue.  The  final  matrix  has  three  spikes,  namely  columns  5,  2  and  7  in  positions 
5,  7  and  8  respectively.  By  interchanpng  colrunns  5  and  7,  we  would  obtsun  only 
two  spikes  but  the  resulting  5x5  leading  submatrix  would  be  structurally  singular. 
The  FIFO  rule  prevents  such  occurences  when  the  whole  matrix  is  non-singular. 


4.6.2  Second  algorithm: 
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Matrix  0 

Iteration  1 

Column  1  is  selected  and  assigned  position  p  =  9  —  2  =  7. 


Rows  1  and  4  axe  selected  and  assigned  positions  7  ^lnd  8. 
Position  8  is  sent  to  the  SP  queue. 
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Iteration  2 


Column  8  is  selected  and  assigned  position  p  =  7  —  2  =  5. 


Rows  7  and  8  are  selected  and  aissigned  positions  5 

Position  6  joins  8  in  the  SP  queue. 
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Matrix  2b 


Iteration  3 

Column  3  is  selected  and  assigned  position  p  =  5  —  1  =  4. 
Row  5  is  selected  and  assigned  position  4. 
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Matrix  3b 


Iteration  4 

Column  2  is  selected  and  assigned  position  p  =  4  —  1  =  3. 
Row  2  is  selected  and  assigned  position  3. 
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Matrix  4b 


Iteration  5 

Column  6  is  selected  and  assigned  position  p  =  3  —  1  =  2. 
Row  6  is  selected  and  assigned  position  2. 
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Matrix  5b 


Iteration  6 

Column  7  is  selected  and  assigned  position  q  =  8. 
Position  8  is  removed  from  the  SP  queue. 
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Iteration  7 

Column  4  is  selected  and  assigned  position  p  =  2  —  1  =  1. 
Row  3  is  selected  and  assigned  position  1. 
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Matrix  7b 

After  7  iterations,  all  rows  have  been  deleted. 

The  remaining  column  5  is 

matched  with  the  position  remaining  in  the  SP  queue.  The  final  matrix  has  only 

two  spikes,  namely  columns  5  and  7,  in  positions  6  and  8  respectively. 

4.6.3 

Third  algorithm: 
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Matrix  0 

Iteration  1 

Column  1 

is  selected  and  assigned  position  p  =  9  —  2 

=  7. 

Rows  1  and  4  are 

selected  and 

assigned  positions  7  and  8. 

Position  8 

is  sent  to  the  SP  queue. 
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Iteration  2 

Row  3  is  selected  and  assigned  position  1. 

Columns  4  ajid  5  are  selected. 

Column  4  is  assigned  position  1  and  column  5  sent  to  the  SI  queue. 
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Iteration  3 

Column  8  is  selected  and  assigned  position  p  =  7  —  2  =  5. 
Rows  7  and  8  are  selected  and  assigned  positions  5  and  6. 
Position  6  is  sent  to  the  SP  queue. 
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Iteration  4 

Row  6  is  selected  and  assigned  position  2. 

Columns  6  and  7  are  selected. 

Column  6  is  assigned  position  2  and  column  7  is  added  to  the  SI  queue. 
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Iteration  5 

Column  3  is  selected  and  assigned  position  p  =  5  —  1  =  4. 
Row  5  is  selected  and  assigned  position  4. 
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Matrix  c5 

Iteration  6 

Row  2  is  selected  and  assigned  position  3. 
Column  2  is  selected  and  assigned  position  3. 
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After  6  iterations,  all  rows  and  all  columns  have  been  deleted.  By  matching 
the  spike-index  queue  {5,7}  with  the  spike-position  queue  {8,6}  in  reverse  order, 
we  obtain  the  same  final  matrix  as  with  the  second  algorithm.  The  two  spikes, 
columns  5  and  7,  are  in  positions  6  and  8  respectively.  In  this  particular  example, 
the  third  algorithm  tjdces  as  few  iterations  as  the  first  one  and  yields  as  few  spikes 
cis  the  second  one. 
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CHAPTER  5.  FACTORIZATION  ALGORITHM 


5.1  Introduction 

The  spike  reduction  algorithms  described  in  Chapter  4  do  not  take  into  con¬ 
sideration  the  numerical  values  of  the  nonzero  coefficients  of  the  matrix  A.  There¬ 
fore,  they  are  unlikely  to  produce  a  matrix  whose  direct  factorization  would  be 
stable.  On  the  other  hand,  the  factorization  with  column  interchanges  described 
in  Chapter  3  would  result  in  too  many  spikes  and  too  much  computation  if  it  were 
applied  to  the  original  matrix  A. 

The  algorithm  described  in  this  chapter  attempts  to  create  a  sparse  and  stable 
factorization  by  combining  the  ideas  described  in  the  previous  three  chapters.  It 
consists  of  three  phases:  prescaling,  preordering  and  factorization. 

The  prescaling  phase  is  a  simple  column  scaling.  The  preordering  phaise  is 
based  on  the  composite  spike  minimization  algorithm  of  Section  4.5.  The  factor¬ 
ization  phase,  which  also  includes  some  reordering  and  rescaling,  is  a  restricted- 
pivoting  version  of  the  factorization  with  column  interchanges  described  in  Section 
3.6. 


5.2  Notation 

We  use  the  definitions  amd  notation  of  Section  4.2.  In  addition,  we  introduce 
a  pivot  tolerance  r  and  the  following  quamtities: 


ROWPivoT(0  =  maoclAjfcl 
coLPivoTO)  =  m^|Ajfc)| 

ROWSCORE(i)  =  f  R0WPIV0T(.)  X  ROWTALLV(.)  if  ROWPIVOT(0  >  T 
\  0  otherwise 

COLSCOREO)  =  /  COLPIVOTO)  X  COLTALLY(j)  if  COLPIVOT(»  >  T 
1 0  otherwise 

Thus,  if  a  row  i  satisfies  ROWSCORE(i)  =  0,  it  contaiins  no  entry  that  can  be  accepted 
ats  a  pivot. 
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5.3  Prescaling  Phase 

To  render  meaningful  any  row  or  column  selection  based  on  the  size  of  the 
coefficients  within  a  column  or  a  row  of  A,  these  coefficients  of  A  must  have  been 
prescaled.  Here,  we  assume  a  columnwise  matrix  representation  and  we  choose  to 
scale  the  columns  of  A.  For  instance,  we  can  divide  each  column  by  a  coefficient 
of  largest  absolute  value  in  that  column.  The  resulting  columns  are  unit  vectors 
for  the  II  .  II norm. 


5.4  Preordering  Phase 

The  preordering  is  aimed  at  reducing  the  number  and  size  of  the  spikes  and, 
hence,  of  the  factors.  It  uses  the  composite  spike  reduction  algorithm  of  Section 

4.5  with  three  modifications. 

The  first  modification  concerns  row  or  column  selection.  The  selected  rows 
and  columns  that  do  not  contain  a  suitable  pivot  (i.e.  an  entry  whose  absolute 
value  is  greater  than  r)  are  rejected.  When  a  column  has  been  selected  and 
several  rows  are  available  to  be  assigned  the  pivot  position,  the  row  containing 
the  entry  of  largest  absolute  value  in  that  column  inside  the  active  submatrix  is 
chosen  as  pivot  row,  unless  that  entry  remains  too  small  to  be  used  as  a  pivot  in 
which  case  the  selected  column  becomes  a  spike.  Similarly,  when  a  row  has  been 
selected  and  several  columns  are  available  to  be  assigned  the  pivot  position,  the 
column  containing  the  entry  of  largest  absolute  value  in  that  row  inside  the  active 
submatrix  is  chosen  as  pivot  colurxm,  unless  that  entry  remains  too  small  to  be 
used  as  a  pivot  in  which  case  the  pivot  position  is  sent  to  the  spike-position  queue 
and  all  the  columns  intersecting  the  selected  row  are  sent  to  the  spike-index  queue. 

The  second  modification  concerns  the  tiebreaking  function  tally  ()  that  is 
replaced  by  the  “multiobjective”  function  scoRE()  to  take  into  account  both  the 
pivot  size  and  the  number  of  nonzero  entries  to  be  deleted. 

Finedly,  the  third  modification  concerns  the  ordering  of  the  spikes.  Spikes  are 
still  given  in  order  of  nonincre2ising  height,  but  their  final  order  is  to  be  determined 
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in  the  factorization  phase.  Therefore,  spike-indices  and  spike-positions  are  only 
added  to,  but  not  removed  from  their  respective  queues. 


5.5  Factorization  Phase 

The  preordering  phase  induces  a  partition  of  the  columns  into  spikes  and 
non-spikes.  The  factorization  phase  preserves  the  ordering  of  the  rows  and  that 
of  the  non-spikes  but  permutes  the  spikes  according  to  the  pivot-maximizing  rule 
introduced  in  Section  3.6. 

For  all  spike-positions  k  in  increasing  order  (i.e.  from  left  to  right),  it  com¬ 
putes  a  vector  Wk  and  the  potential  pivots  associated  with  the  available  spikes, 
selects  a  pivot  of  largest  absolute  value  and  matches  the  corresponding  spike-index 
with  the  spike-position  k.  If  no  suitable  pivot  can  be  obtained  from  the  spike-index 
queue,  the  matrix  is  deemed  singular  and  the  algorithm  stops. 

In  addition,  each  time  a  spike  has  been  assigned  a  position,  it  may  be  rescaled 
so  that  the  resulting  pivot  takes  the  value  1.  This  step  is  not  recommended  for 
dense  matrices  (cf.  Section  3.7). 


5.6  Factorization  Algorithm 

All  the  procedvures  mentioned  above,  including  the  column  rescaling,  are  con¬ 
tained  in  the  following  algorithm: 


Normalize  the  columns  of  A. 

Let  fc  =  m  +  1. 

Let  pfc  =  m  +  1. 

Let  k'  =  0. 

Repeat 

If  not  all  columns  are  deleted  then 
Let  k  =  k  —  1; 
bottomright(  Jfc) . 

If  not  all  rows  axe  deleted  then 
Let  k'  =  k'  +  1; 
topleft(A:'). 
until  all  rows  are  deleted. 

Add  imdeleted  columns  to  the  spike-index  queue. 

MATCHQUEUES. 

The  procedure  BOTTOMRiGHT(fc)  becomes: 

•  Select  from  the  active  submatrix  a  column  jk  minimizing  colcount(»;  breeik 
eventual  tie  by  maximizing  coLSCORE(j). 

•  Update  the  rightmost  non-spike  position  pk  =  Pk+i  —  C0LC0UNT(». 

•  If  COLSCOREO*)  ^  0,  then 

assign  column  jk  position  pk; 

select  a  pivot  row  from  the  rows  intersecting  column  jk‘, 
assign  the  pivot  row  position  pk; 

assign  the  other  rows  intersecting  column  jk  the  positions  in  the  range 
p*-|-l,...,pjk+i  -1; 

swid  the  positions  p*  -|- 1 . .  .pfc+i  —  1  to  the  spike-position  queue; 
rescale  the  pivot  column  so  that  the  resulting  pivot  equals  1. 

•  If  coLSCORE(j»)  =  0,  then 

add  column  jk  to  the  spike-index  queue. 

•  Delete  column  jk  and  the  rows  intersecting  column  jk- 
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The  procedure  TOPLEFT(fc)  becomes: 

•  Select  from  the  active  submatrix  a  row  minimizing  ROWcouNT(i)  ;  break 
eventual  tie  by  maximizing  ROWSCORE(i)- 

•  Assign  row  ik  position  k. 

•  If  ROWScoRE(u)  ^  0,  then 

select  a  pivot  column  jk  from  the  columns  intersecting  row  ik] 
aissign  column  jk  position  k] 

add  the  other  columns  intersecting  row  ik  to  the  spike- index  queue. 

•  If  R0WSC0RE(i*)  =  0  then 

add  position  k  to  spike-position  queue; 

add  all  columns  intersecting  row  ik  to  spike-index  queue. 

•  Delete  row  ik  and  the  columns  intersecting  row  i*. 

The  procedure  matchqueues  stands  for: 

•  While  spike-position  queue  0, 

remove  leftmost  position  k  from  spike-position  queue; 
compute 

compute  the  potential  pivots  associated  with  the  columns  of  the  spike- 
index  queue; 

select  pivot  ^k  and  pivot-maximizing  column  /; 

if  l^fcl  <  e,  then  STOP  (A  is  singular); 

assign  /  position  k  and  remove  \  from  spike-index  queue; 

rescale  the  pivot  column  so  that  the  resulting  pivot  equals  1. 


5.7  Example: 

In  this  Section,  we  apply  the  sparse  factorization  algorithm  described  above 
to  the  8x8  matrix  introduced  in  chapter  4.  The  pivot  tolerance  r  is  set  to  0.50. 

I 

The  numerical  values  of  the  matrix  and  of  its  factors  are  given  with  a  precision  of 
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10“^.  Phases  1,  2  and  3  correspond  to  prescaling,  preordering  and  factorization 
with  restricted  pivoting  respectively.  SI  queue  and  SP  queue  denote  the  spike- 
index  queue  and  the  spike-position  queue  respectively.  The  columns  belonging  to 
the  SI  queue  axe  printed  in  italic. 
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Phase  2  bottomrightI 

Column  1  is  selected  and  assigned  position  p  =  9  —  2  =  7. 
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Rows  1  (pivot  row)  and  4  are  selected  and  assigned  positions  7  and  8. 
Position  8  is  sent  to  the  SP  queue. 
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Phase  2  topleftI 

Row  3  is  selected  and  assigned  position  1. 
Columns  4  and  5  are  selected. 

Column  4  (pivot  column)  is  assigned  position  1. 
Coliunn  5  is  sent  to  the  SI  queue. 

Colmnn  4  is  rescaled  so  that  the  pivot  equals  1. 
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Column  3  is  selected  and  assigned  position  p  =  7  —  2  =  5. 

Rows  5  (pivot  row)  and  8  are  selected  and  assigned  positions  5  and  6. 
Position  6  is  sent  to  the  SP  queue. 
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Phase  2  topleft2 

Row  6  is  selected  and  assigned  position  2. 
Columns  6  and  7  are  selected. 

Column  6  (pivot  column)  is  eissigned  position  2. 
Column  7  is  sent  to  the  SI  queue. 

Column  6  is  rescaled  so  that  the  pivot  equeds  1. 
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Phase  2  bottomrightS 

Column  8  is  selected  and  assigned  position  p  =  5  —  1  =  4. 
Row  7  (pivot  row)  is  selected  and  assigned  position  5. 
Column  8  is  rescaled  so  that  the  pivot  equals  1. 
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Phase  2  topleft3 

Row  2  is  selected  and  rejected  because  |0.33|  <  r. 
Position  3  is  sent  to  the  SP  queue. 

Column  2  is  sent  to  the  SI  queue. 


Phase  3  Computation  of  ^3 

We  remove  the  leftmost  position,  namely  3,  from  the  spike-position  queue.  Then 
we  have  to  solve  the  system  u^sAi  =  uj’,  or  rather 


/1.00  \ 

^3  1.00  I  =  (0  0  1). 

V  1.60  1.00/ 

The  solution,  7^3  =  ( 0.00  —1.60  1.00 ),  can  be  used  to  compute  the  pivots  that 

would  result  from  assigning  columns  2,  7  and  5  position  3 

/  -0.50  \ 

(0.00  -1.60  1.00)  1  -0.17  0.25  |  =  (0.33  0.60  -0.40). 

\0.33  0.33  / 
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Selecting  column  7  yields  the  pivot  of  largest  absolute  veilue.  Therefore,  we  remove 
column  7  from  the  spike-index  queue  and  assigned  it  position  3.  Then  we  rescale 
column  7  by  0.60  so  that  the  resulting  pivot  equals  1. 
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VVe  remove  the  leftmost  position,  6,  from  the  spike-position  queue.  Then  we  have 
to  solve  the  system  cjgEaAi  =  uj ,  or  rather 
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1.00),  can 

be  used  to 

compute  the  pivots  that  would  result  from  assigning  columns  5  and  2  position  6: 
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Column  5  provides  the  pivot  of  largest  absolute  vsJue.  Therefore,  we  remove 
column  5  from  spike-index  queue  and  assign  it  position  6.  Then  we  rescale  it  by 
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1.59  so  that  the  resulting  pivot  equ2ils  1. 
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-7.00 

-1  2.7 

8.00 

9  00 

omega 

-0.77 

1.83 

-1.33 

0.75 

0.29 

1.00 

beta 

1  59 

0  24 

Matrix  8d 

Phase  3  Computation  of  (Ta 

We  match  the  remaining  spike-index  2  and  spike-position  8.  However,  we  must 


still  solve  oJgEgEgAi  =  uj  or  equivalently 

/  1.00 

-0.31 

1.00  -0.28 

0.16 

1 

1.60  0.56 

0.83 

0.40  1.67 

1.00 

0.63 

0.50 

1.00 

-0.75  - 

-0.29 

0.17 

1.75  - 

-0.71 

1.00 

1 

\-1.50 

-0.50 

-0.38 

1.00  y 

=  (0  0 

0  0  0 

0  0 

!)• 

The  solution  is  yg  =  (0.97  0.73  -0.53  0.30  0.44  0.61  0.38  1.00) 
and  the  pivot  resulting  from  placing  column  2  in  position  8  is  given  by 


(0.97  0.73  -0.53  0.30  0.44  0.61  0.38  1.00) 


(  \ 

0.33 

0.44 

0.56 

1.00 

\0.56/ 


=  1.29 
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Therefore,  we  rescale  column  2  by  1.29  so  that  the  pivot  equals  1. 


4  6  7 


3 

1  00 

6 

1.00 

-0  28 

2 

1.60 

0.56 

7 

0  83 

0.40 

1  67 

5 

0  50 

8 

1 

0  17 

4 

-  1.50 

scaie 

6.00 

5.00 

3  60 

omega 

0.97 

0  73 

-0.53 

beta 

8 

3 

5 

1 

2 

-0  3  1 

0.16 

0.2C 

1.00 

1.00 

0.63 

0.35 

-0  75 

-0.29 

0.43 

1.75 

-0.71 

1.00 

0.78 

-0.50 

-0.38 

0.43 

-4.00 

-7  00 

-12.7 

8.00 

1  1  58 

0.30 

0.44 

0  6  1 

0.38 

1 .00 

1.29 

Matrix  9d 


The  final  matrix  has  three  spikes,  namely  columns  7,  5  and  2,  in  positions  3, 
6  and  8  respectively.  Its  factorization  is  given  by  the  three  vectors 
^3  =(0.00  -1.60  1.00) 

^6  =  (-0.77  1.83  -1.33  0.75  0.29  1.00) 

^8  =  (0.97  0.73  -0.53  0.30  0.44  0.61  0.38  1.00). 


5.8  Block  Triangular  Reduction 

Dulmage  and  Mendelsohn  (1963)  have  indicated  a  procedure  to  permute  the 
rows  and  the  columns  of  a  matrix  so  that  the  resulting  matrix  is  lower  triangular 
by  blocks. 

First,  find  a  transversal,  i.e.  a  permutation  of  the  columns  with  matrix  rep¬ 
resentation  Q  such  that  the  diagonal  of  the  resulting  matrix  AQ  has  no  zero 
entry. 

Then  identify  the  resulting  matrix  AQ  with  its  canonically  associated  directed 
graph  (the  nodes  are  the  indices  {!,...  ,m}  and  the  au’cs  axe  the  pairs  {i,j)  such 
that  Ax,j  ^  0). 
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Finally  determine  the  strongly  connected  components  of  the  graph  (cf.  Tarjan, 
1972),  regroup  the  nodes  of  the  same  components  into  supemodes  and  consider 
the  resulting  collapsed  graph:  since  it  does  not  contain  ciny  cycle,  the  supernodes 
can  be  reordered  so  that  if  arc  /  -+  J  exists  then  I  >  J .  List  the  nodes  from  the 
first  supernode,  then  those  from  the  second  one,  gind  so  on  that  list  defines  a  row 
permutation  matrix  P.  The  matrix  PAQP^  is  lower  triangular  by  blocks. 

Once  the  blocks  have  been  identified,  they  can  be  reordered  so  as  to  minimize 
the  number  of  spikes  within  the  block  and  systems  can  be  solved  block  by  block. 
This  reduces  the  size  of  the  blocks  to  be  reordered  and  the  length  of  the  auxiliary 
vectors  or  (Tk- 
Example; 


On  the  10  X  10  matrix  shown  above,  factorizing  the  diagonal  blocks  instead 
of  the  whole  matrix  reduces  the  number  of  nonzero  components  for  the  auxiliary 
vectors  ufk  or  (Tk  from  25  to  17. 
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CHAPTER  6.  APPLICATION  TO  THE  SIMPLEX  ALGORITHM 


6.1  Introduction 

In  order  to  test  our  factorization  method,  we  have  implemented  it  within  the 
framework  of  the  simplex  algorithm  for  lineair  programming  (Dantzig,  1963).  As 
a  reference  code,  we  have  chosen  the  FORTRAN  optimization  code  MINOS  5.3 
(Murtagh  and  Saunders,  1987)  whose  modularity,  robustness  and  performance  on 
leuge-scEile  problems  malce  it  an  ideal  benchmark. 

In  MINOS  5.3,  the  factorization  of  the  basis,  the  update  of  the  basis  and  the 
solution  of  linear  systems  involving  the  basis  are  cauried  out  by  a  set  of  FORTRAN 
subroutines  constituting  the  file  MI25BFAC.  We  have  written  a  file  of  FORTRAN 
subroutines,  called  MI26BFAC,  designed  to  perform  the  same  taisks  using  our 
method.  We  have  run  MINOS  5.3  alternatively  with  the  original  file  MI25BFAC 
and  with  the  new  file  MI26BFAC  on  different  test-problems  under  different  options. 


6.2  The  File  MI26BFAC 

The  file  MI26BFAC  is  designed  to  perform  the  same  tasks  as  MI25BFAC 
when  used  within  MINOS  5.3  to  solve  lineau:  programming  problems.  The  major 
implementation  differences  are  the  following.  First,  our  factorization  is  inherently 
different  from  the  LU  factorization.  Second,  we  compute  a  block-triangularization 
of  the  basis  and  then  perform  the  factorization  only  on  the  diagonal  blocks.  Third, 
we  use  a  conventional  product-form  update  (with  an  in-core  “7/-file”). 

The  principal  subroutines  of  MI26BFAC,  written  in  FORTRAN??,  are  repre¬ 
sented  hierarchically  in  Chart  6.2.  Their  names  and  their  roles  are  listed  below: 
M2BFAC:  Performs  the  factorization  of  the  beisic  matrix  B  by  calling  N2BELM, 
M2BMAP,  M2BS0L  and  M2SING. 

M2BELM:  Extracts  a  sequence  of  triplets  representing  the  matrix  B  from  the  data. 
POINTR:  Converts  the  representation  of  B  from  a  sequence  of  triplets  to  sparse 
colunm  format  with  an  array  of  coefficient  values,  an  array  of  row  indices 
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and  2m  array  of  column  pointers. 

M2BMAP:  Allocates  storage  for  most  of  the  arrays  used  in  the  factorization  and 
updating  of  B. 

M2BS0L:  Performs  one  of  the  following  actions  according  to  the  value  of  the  pa¬ 
rameter  “mode”; 

factorize  B  by  calling  POIHTR,  TRIAHG  and  TRTBLK; 
solve  BHi  . .  .HpPx  =  b  by  calling  SOLVEM,  PERMCP  and  UPDATE; 
solve  ttBHi  . . .  Hp  =  7P-*  by  calling  SOLVER,  PERMCP  and  UPDATE; 
add  an  element  H,  to  the  rj-file  by  calling  ADDETA. 

H2SIHG:  Replaces  one  coltimn  of  B  by  an  appropriate  slack  column  when  the 
original  matrix  results  in  too  small  a  pivot  during  factorization.  It  may 
be  called  several  times  in  a  row. 

TRIAHG:  Identifies  the  lower  block  triangular  representation  of  B  by  calling  TRAHSV 
and  GETBLK. 

TRTBLK;  Extracts,  prescales,  preorders  and  factorizes  each  diagonal  block  of  B  by 
calling  EXTRCT,  RESCAL,  PREORD  and  FACTOR. 

SOLVEM:  Solves  Bx  =  b  with  the  help  of  SOLVEA. 

SOLVEH:  Solves  irB  =  7  with  the  help  of  SOLVEB. 

PERMCP:  Permutes  a  vector  according  to  P  or  P"*. 

UPDATE:  Performs  rank-one  updates  according  to  the  T7-file. 

ADDETA:  Adds  a  coliunn-vector  characterizing  a  factor  Hj  to  the  T/-file. 

TRAHSV:  Computes  a  transversal  of  B. 

GETBLK:  Computes  the  diagonal  blocks  of  B. 

EXTRCT:  1.  ^c^ntifies  the  elements  from  a  diagonal  block  of  B  with  the  help  of  SELECT. 

SELECT:  Scans  the  coliunns  that  support  a  diagonal  block  to  identify  the  entries 
that  belong  to  that  block. 

RESCAL:  Normalizes  the  columns  of  each  block. 

PREORD:  Applies  a  spike  reduction  algorithm  to  preorder  the  blocks  by  calling 
SETLST,  COLCAH,  BOTRIG,  ROWCAH,  TOPLEF,  PRUHE  and  GATHER. 
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FACTOR:  Factorizes  the  blocks  by  calling  SOLVES  and  ADDSPK.  Sorts  the  column 
representation  of  B  according  to  the  final  row  order. 

SETLST:  Sets  up  a  row-wise  linked  list  representing  a  diagonal  block  and  some 
auxiliary  variables  for  preordering. 

COLCAN:  Finds  colunrn  candidates  to  be  used  by  BOTRIG. 

BOTRIG:  Performs  a  BottomRight  iteration. 

ROWCAN;  Finds  row  candidates  to  be  used  by  TOPLEF. 

TOPLEF:  Performs  a  TopLeft  iteration. 

PRUNE  :  Performs  the  deletion  of  entries  in  row-wise  linked  lists. 

GATHER:  Gathers  and  orders  all  spikes  in  order  of  non- increasing  height. 

SOLVEA:  Solves  B'x'  =  b'  where  B'  is  a  sub-block  of  B. 

SOLVEB:  Solves  ir'B'  =  7'  where  B'  is  a  sub-block  of  B. 

ADDSPK:  Adds  one  auxiliary  vector  tTk  to  the  factor  list. 

6.3  Test  Implementation 

We  have  run  MINOS  5.3  on  a  SUN  3/50  workstation  using  a  16  Mhz  Motorola 
68020  CPU  under  the  SunOS  3.5  operating  system  and  with  the  Berkeley  f77 
compiler. 

We  have  chosen  a  set  of  53  test- problems  studied  by  Lustig  (1987)  and  made 
publicly  available  on  aetlib  (Dongarra  and  Grosse,  1987)  by  Gay  (1985).  Although 
the  original  MINOS  code  could  run  adl  these  test-problems  without  modification, 
our  implementation  required  too  much  storage.  Therefore,  we  have  increased 
the  size  of  the  *irray  Z  from  100000  to  120000  for  SHIP12L  and  to  150000  for 
PILOTJA  and  decided  to  forego  the  last  two  problems  80BAU3B  and  PILOT 
because  it  was  appeurent  from  running  the  other  test-problems  that  our  method 
would  be  unpractical  and  widely  outperformed  by  MINOS. 

With  each  of  the  51  other  test-problems,  we  have  made  two  experiments,  each 
based  on  a  different  MINOS  option  file. 
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6.4  First  Experiment 


The  first  experiment  is  the  straightforward  solve  of  each  lineeir  program  from 
scratch  in  order  to  compare  overall  performance.  The  option  file  is  the  following: 


ROWS 

1500 

COLUMNS 

5500 

ELEMENTS 

22000 

MPS  FILE 

10 

NEW  BASIS  FILE 

12 

SAVE  FREQUENCY 

10000 

ITERATION  LIMIT 

20000 

PRINT  LEVEL 

0 

SOLUTION 

NO 

SCALE  OPTION 

2 

PARTIAL  PRICE 

10 

FACTOR  FREQUENCY 

25 

LU  FACTOR  TOLERANCE 

100.0 

LU  UPDATE  TOLERANCE 

10.0 

In  Table  6.4.1,  we  have  listed  for  each  test-problem  the  number  of  iterations 
/„  in  phase  1,  the  total  number  of  iterations  and  the  execution  time  under 
both  versions  (u  =  0  for  MI25BFAC  and  u  =  1  for  MI26BFAC),  as  well  zis  the 
ratio  of  execution  times  ti/to-  In  19  cases,  our  implementation  runs  faster  than 
MINOS,  which  is  hampered  by  the  factorization  frequency  of  25. 

In  Table  6.4.2,  we  have  listed  for  each  test-problem  the  optimal  objective 
value  Zv  under  both  versions  (u  =  0  for  MI25BFAC  and  u  =  1  for  MI26BFAC), 
as  well  as  their  relative  difference.  In  44  instances,  the  optimal  objective  vzdues 
agree  up  to  11  significant  digits.  The  largest  discrepancy  occurs  with  PILOTWE 
where  the  agreement  is  still  of  6  signific2mt  digits.  This  seems  to  indicate  that  our 
method  achieves  a  satisfactory  numerical  stability. 


71 


kO 

_ t1 

to 

rti/to 

ARRO 


ADUTTLE 


SC205 


SCAGR7 


SHARE2B 


RECIPE 


VTPBASE 


SHARE1B 


BORE3D 


SCORPION 


CAPRI 


SCAGR25 


SCTAP1 


BRANDY 


ISRAEL 


ETAMACflO 


SCFXM1 


GRQW7 


BANDM 


E226 


STANOATA 


SCS01 


GFRDPNC 


BEAOOMTI 


STAIR 


SCRS8 


SEBA 


SHELL 


PILOT4 


SCFXM2 


SCS06 


GROW15 


SHIP04S 


FFFFF800 


QYces 


SCFXM3 


SCTAP2 


GROW22 


SHIP04L 


PILOTWE 


SIERRA 


SHIP08S 


SCTAP3 


SHIP12S 


25FV47 


SCS08 


NESM 


2046 


767 


1176 


1  3 


492 


452 


1  7 


425 


39 


2304 
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1251 


1  10 


33 


4  8 


261 


160 


102 


251 


353 


287 


423 


250 


649 


389 


232 


498 


480 


114 


368 


6151 


1  04 


449 


692 


399 


274 


1599 


819 


1139 


485 


158 


1074 


700 


1349 


766 


736 


275 


4474 


1083 


262 


926 


434 


7821 


3102 


3270 


110 


33 


4  8 


261 


160 


102 


251 


3761 


287 


423 


250 


574 


389 


232 


498 


467 


1  1  1 


370 


637 


8  7 


449 


645 


399 


274 


1467 


819 


1139 


485 


159 


1002 


705 


1391 


744 


662 


276 


40401 


1071 


262 


883 


434 


7828 


3914 


3399 


1.94 


8.00 


35.58 


11.50 


10.52 


5.88 


10.10 


27.86 


27.38 


23.40 


41 .32 


81 .88 


46.78 


80.42 


49.28 


129.30 


70.54 


54.86 


114.10 


81.60 


31.001 


143.021 


25.62 


385.68 


176.42 


1 12.00 


73.26 


1278.24 


259.50 


172.96 


226.32 


47.12 


292.22 


333.02 


577.54 


319.02 


571.92 


77.60 


7900.80 


438.18 


106.68 


488.96 


196.14 


10704.38 


1383.04 


1357.34 


1.84 


7.18 


17.84 


10.14 


9.48 


6.30 


9.72 


23.56 


25.46 


24.88 


37.72 


86.00 


47.72 


65.66 


37.18 


109.30 


68.98 


49.82 


101 .06 


67.76 


31.10 


37.68 


7.26 


22.82 


197.58 


161 .68 


106.14 


82.02 


550.82 


270.40 


148.24 


196.24 


54.42 


284.08 


354.66 


660.18 


346.84 


369.46 


88.86 


1938.60 


511.84 


130.24 


546.74 


273.56 


4524.18 


1245.60 


1189.00 
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TABLE  6.4.2 

Zl 

-  - 

AFRO 

-4.64753142869+02 

ADLITTLE 

2.25494963169+05 

SC205 

-5.22020612129+01 

SCAGfl? 

-2.33138925489+06 

SHARE2B 

-4.15732240749+02 

RECIPE 

-2.66616000009+02 

VTPBASE 

1.29831462469+05 

SHARE1B 

-7.65893185799+04 

B0RE30 

1.37308039429+03 

SCORPION 

1.87812482279+03 

CAPRI 

2.69001291389+03 

SCAGR25 

-1.47534330619+07 

SCTAP1 

1.41225000009+03 

BRANDY 

1.51850989659+03 

ISRAEL 

-8.96644821869+05 

ETAMACRO 

-7.55715218199+02 

SCFXM1 

1.84167590289+04 

GROW? 

-4.778781 18159+07 

8ANDM 

-1.58628018459+02 

E226 

-1.87519290669+01 

STANDATA 

1.25769950009+03 

SCSD1 

8.66666667439+00 

GFFDPNC 

6.90223599959+06 

BEAOOMD 

3.35924858079+04 

STAIR 

•2.51266951199+02 

SCRS8 

9.04299986199+02 

SEBA 

1.57116000009+04 

SHELL 

1.20882534609+09 

PIL0T4 

-2.581 13926419+03 

SCFXM2 

3.66602615659+04 

SCS06 

5.05000000789+01 

GROW15 

•1.06870941299+08 

SHIP04S 

1.79871470049+06 

FFFFF800 

5.55679612199+05 

GANGES 

-1.09586347709+05 

SCFXM3 

5.49012545509+04 

SCTAP2 

1.72480714299+03 

GROW22 

1.60834336489+08 

SHIP04L 

1.79332453809+06 

PILOTWE 

-2.7200991 1969+06 

S€FIRA 

1.53943621849+07 

SHIP08S 

1.92009821059+06 

SCTAP3 

1 .42400000009+03 

SHIP12S 

1.48923613449+06 

25FV47 

5.50184677919+03 

scsoe 

9.04999999939+02 

NESM 

1.40760733249+07 

CZPROB 

2.16519669699+06 

PILOTJA 

•6.1  1311529489+03 

SHIP08L 

1.90905521149+06 

SHIP12L 

1.47018791939+06 

•4.6475314286e-t-02 


2.254949631  6e-»-05 


-5.2202061  21  2e-t-01 


•2.331  3892548e-i-06| 


-4.15732240749+02 


-2.66616000009+02 


1.29831462469+05 


-7.65893185799+04 


1 .37308039429+03 


1.87812482279+03 


2.69001291389+03 


-1.47534330619+07 


1.41225000009+03 


1.51850989659+03 


-8.96644821869+05 


-7.55715218319+02 


1 .84167590289+04 


-4.77878118159+07 


-1.58628018459+02 


-1.87519290669+01 


1.25769950009+03 


8.66666667439+00 


6.90223599959+06 


3.35924858079+04 


-2.51266951199+02 


9.04299986199+02 


1.57116000009+04 


1.20882534609+09 


-2.581 13926419+03 


3.66602615659+04 


5.05000000789+01 


-1.06870941299+08 


1.79871470049+06 


5.55679591039+05 


•1.09586352259+05 


5.49012545509+04 


1.72480714299+03 


1.60834336489+08 


1.79332453809+06 


-2.72010345259+06 


1.53943621849+07 


1.92009821059+06 


1.42400000009+03 


1.48923613449+06 


5.50184588839+03 


9.04999999939+02 


1.40760751289+07 


2.18519669699+06 


-6.1 131 1577759+03 


1.90905521 149+06 


1.47018791939+06 


zl-zO/zO 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


-1 .599-1  0 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


3.819-08 


-4.159-08 


0.009+00 


0.009+00 


0.009+00 


0.009+00 


-1.599-06 


0.009+00 


0.009+00 


0.009+00 


1.629-07 


0.009+00 


-1.289-07 


0.009+00 


-7.909-08 


0.009+00 


0.009+00 
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6.5  Second  Experiment 

The  second  experiment  is  designed  to  focus  on  the  factorization  size.  As  initial 
bcisis,  we  use  the  optimal  basis  obtained  during  the  first  experiment  under  MINOS 
with  MI25BFAC.  The  objective  function  is  now  maximized  instead  of  minimized 
and  the  iteration  limit  is  set  to  25.  The  option  file  is  the  following: 


ROWS 

1500 

COLUMNS 

5500 

ELEMENTS 

22000 

MAXIMIZE 

MPS  FILE 

10 

OLD  BASIS  FILE 

12 

NEW  BASIS  FILE 

0 

SAVE  FREQUENCY 

10000 

ITERATION  LIMIT 

25 

PRINT  LEVEL 

1 

SOLUTION 

NO 

SCALE  OPTION 

2 

PARTIAL  PRICE 

10 

FACTOR  FREQUENCY 

100 

LU  FACTOR  TOLERANCE 

100.0 

LU  UPDATE  TOLERANCE 

10.0 

Not  surprisingly,  some  of  the  test-runs  end  with  an  unbounded  solution.  In 
four  other  instances,  MINOS  lists  some  variables  as  apt  to  increase  indefinitely. 
These  occurences  are  indicated  respectively  by  a  U  and  an  I  in  the  first  column  of 
Table  6.5.1. 

Table  6.5.1  also  contains  the  number  of  iterations  k  up  to  which  both  methods 
yield  the  same  basic  solution,  and  the  total  number  of  iterations  k' .  For  each  test- 
problem,  that  number,  usually  the  iteration  limit  of  25,  turns  out  to  be  the  same 
under  both  versions. 

Finally,  Table  6.5.1  provides  the  size  cl  of  the  largest  diagonal  block,  the 
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m\ 


25 

wm 

■ElBE] 

8 

2 

1 

_ 5 

0.58 


0.28 


0.12 


0.08 


0.12 


0.08 


12.36 


0.34 


0.14 


0.26 


9.42 


0.58 


0.18 


2.10 


0.16 


0.32 


0.98 


0.86 


0.38 


4.08 


0.18 


15.70 


0.22 


0.24 


0.40 


0.36 


6.82 


1.14 


1.08 


0.30 


25.86 


0.24 


0.30 


1.00 

1.70 

1.98 

0.86 

0.43 

4.32 

4.22 

1.02 

1.26 

11.06 

7.38 

1.50 

0.22 

5.72 

5.50 

1.04 

0.42 

5.70 

5.34 

1 .07 

0.1 1 

5.40 

6.32 

0.85 

0.28 

7.88 

7.84 

1.01 

0.29 

7.62 

7.50 

1.02 

0.41 

1 1 .52 
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1.02 

0.32 
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0.95 

0.34 

13.22 
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1.11 

0.30 

16.14 
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1 .05 
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12.86 

1.03 

0.67 
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1.11 

0.37 

8.94 

9.24 
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0.30 
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17.50 

1.00 
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15.42 

1 .01 

0.96 

18.04 
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1.17 

0.48 
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1.06 

0.39 
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1.08 
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0.99 

0.36 

15.00 
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26.40 
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0.21 
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12.06 

1 .00 

2.75 
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1.73 
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1.00 

0.1 1 
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1.02 

0.29 
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32.10 
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25.18 
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0.54 
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1.01 
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0.96 
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0.96 
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42.14 
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54.72 

1.33 

0.75 

52.32 

51.30 

1.02 
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4.34 
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3424 

551 
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fart orizat ion  times  /,,  under  both  versions  (u  =  0  for  MI25BFAC  and  u  =  1  for 
Mr2GBFAC),  as  well  as  their  ratio  /;/ /q,  and  the  total  running  time  under  both 
v(>rsions  { r  =  0  for  MI25BFAC  and  =  1  for  MI26BFAC),  as  well  as  their  ratio 
t]/to-  Our  method  outperforms  MINOS  in  44  of  the  factorization  times  but  only 
in  21  of  the  total  execution  times. 

In  addition  to  the  auxiliary  vectors  (or  tT*,-),  our  factorization  requires  the 
storage  <tf  some  elements  of  the  original  basis.  The  only  elements  of  the  basis  that 
iK'od  not  be  stored  rue  located  in  the  lower  triangular  part  of  the  diagonal  blocks. 
Mor<'  prt'cisely,  the  first  k  elements  of  the  rows  of  index  k  where  column  k  is  a 
spik<'  are  not  needed  (f)nce  the  vector  at*  has  been  computed).  We  denote  by  a 
the  numb(u'  of  nonzero  elements  in  the  original  basis,  by  a  the  number  of  these 
nonzero  elements  that  can  he  viewed  as  part  of  the  factorization,  by  w  the  number 
of  nonzero  ('lements  in  the  auxiliary  vectors  utt,  by  x  the  number  of  elements  stored 
in  the  jy-file.  Thus  the  amount  of  storage  required  by  our  method  is  S]  =  k+w  for 
the  original  basis  and  s)  =  a+u)  +  x  for  the  basis. 

Finally,  we  denote  by  l,u  and  I' ,u'  the  sizes  of  the  LU  factors  computed  by 
MINOS  for  the  first  and  basis  respectively.  The  amount  of  storage  required 
by  MINOS  is  So  =  I  +  u  for  the  original  basis  and  s'q  =  I'  +  u'  for  the  basis. 

In  Table  6.5.2,  we  have  listed  the  values  of  a  and  w,  the  sizes  S]  and  s\  of 
the  factorizations  under  MI26BFAC,  the  values  I  and  u,  the  sizes  Sq  and  .Sq  of  the 
factorizations  under  Mr25BFAC,  as  well  as  the  ratios  S]/.:o  and  -Sj/sq.  In  test- 
Itroblerns  where  unboundedness  has  been  detected,  some  factorizations  have  not 
tx  en  carried  out,  nmdering  .some  of  the  above  values  unavailable.  These  instances 
tire  indic;it<’fl  by  a  *. 

In  terms  of  sparsity,  our  method  performs  as  well  as  or  slightly  better  than 
MIXOS  in  32  instances  of  direct  factorizations.  However  it  performs  uniformly 
worse  after  20  column  updates. 
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6.6  Conclusion 

The  testing  experiment  of  Section  6.4  indicates  that  our  factorization  method 
is  numerically  stable  in  practically  all  instances. 

The  testing  experiment  of  Section  6.5  demonstrates  that  our  factorization 
method  is  more  efficient  when  the  diagonal  blocks  computed  during  the  block- 
trianguiarization  are  small,  say  of  order  50  or  less.  This  can  be  attributed  to  the 
following  points: 

•  The  identification  of  the  diagonal  blocks  resulting  from  the  block-triangular- 
ization  is  efficient. 

•  Only  the  diagonal  blocks  are  factored. 

•  A  sorting  of  the  rows  that  speeds  up  the  reordering  of  each  block  is  imple¬ 
mented  along  with  the  factorization. 

On  the  other  hand,  the  same  experiment  shows  that,  when  the  block-triangul- 
arization  yields  a  large  diagonal  block,  the  size  of  our  factorization  becomes  pro¬ 
hibitively  large  and  renders  the  whole  method  inefficient.  In  practically  all  in¬ 
stances  the  numerical  stability  remains  satisfactory. 

Unfortunately,  there  does  not  seem  to  exist  an  efficient  updating  algorithm  for 
the  block-triangularization  of  a  matrix.  Because  our  method  only  factorizes  the 
diagonal  blocks,  this  makes  a  conventional  product-form  updating  almost  manda¬ 
tory.  This  method  of  updating  appears  less  efficient  than  the  Bartels-Golub  up¬ 
dating  of  the  LU  factors  implemented  in  MINOS  5.3. 

In  spite  of  these  drawbacks,  we  hope  that  future  computer  architectures  and 
numerical  software  will  suggest  more  applications  for  all  or  part  of  our  factorization 
method. 
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As  an  alternative  to  the  LU  matrix  factorization,  we  consider  a  factoriza¬ 
tion  that  uses  the  lower  trizingular  part  of  the  originad  matrix  ais  one  factor  and 
computes  the  other  factors  as  a  product  of  rank -one  update  matrices. 

Under  some  non-singularity  assumptions,  an  m  xm  matrix  A  can  be  factorized 
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one  update  matrix  of  the  form  I  -I-  VfcWjt  with  v*  a  column  vector  and  a  row 
vector.  The  vector  v*  is  the  fc''*  column  of  A  —  Aj .  If  v*  =  0,  then  E^  =  I  may  be 
omitted  from  the  factorization.  Otherwise,  the  row  vector  Ufk  must  be  computed. 

After  reviewing  and  improving  the  time  complexity,  the  requirements,  the 
stability  and  the  efficiency  of  this  method,  we  derive  a  stable  factorization  algo¬ 
rithm  which  we  implement  in  FORTRAN??  within  the  fraunework  of  the  simplex 
algorithm  for  linear  programming. 

A  comparison  of  our  numerical  results  with  those  obtained  through  the  code 
MINOS  5.3  indicate  that  our  method  may  be  more  efficient  than  am  ordinary 
LU  decomposition  for  some  matrices  whose  order  ramges  between  28  and  1481, 
especially  when  these  matrices  are  admost  triangular. 
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